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SECTION  I 
INTRODUCTION 


HULL  is  a  system  of  computer  programs  which  solve  the  finite-difference 
analogs  of  the  partial  differential  equations  descriptive  of  a  non-conducting, 
inviscid,  elastic/plastic  continuous  medium  subjected  to  dynamic  loading. 

HULL  is  capable  of  running  cylindrical  or  cartesian  coordinate  calcula¬ 
tions  in  two  spatial  dimensions  or  cartesian  coordinate  calculations  in  three 
dimensions.  HULL  is  a  multimaterial  Eulerian  code  which  employs  a  diffusion 
limiter  to  preseve  material  interfaces.  Up  to  ten  materials  can  bo  handled  in 
a  single  calculation.  A  mixture  of  any  of  these  materials  is  permitted  in  any 
one  cell=  Solid  material  properties  are  modeled  by  tbe  Mie-Gruneisen  equation 
of  state  (REF  1).  Plasticity  is  modeled  by  the  Von  Miees  yield  criteria  which 
can  be  modified  by  plastic  work  hardening  and  thermal  softening.  Slip  sur¬ 
faces  are  automatically  handled  at  material  and  failure  interfaces.  Failure 
can  be  determined  by  simple  stress,  strain  or  tri-axial  criteria.  Failure  is 
modeled  by  inserting  a  void  or  low  pressure  gas  in  the  failed  zone  and  allow¬ 
ing  the  tensile  stresses  to  relax.  Recompression  of  the  failed  zone  is  not 
inhibited. 

Large  computational  meshes  can  be  handled  by  HULL  since  only  a  small 
portion  of  the  mesh  need  be  in  central  memory  as  the  calculation  proceeds. 
Time  snap-shots  at  arbitrary  intervals  can  be  retained  on  permanent  disk  or 
tape  files  for  subsequent  analysis  or  problem  restart.  The  HULL  system  has 
several  service  routines  tor  keeping  track  of  this  data  storage  and  for  pro¬ 
viding  graphical  presentati  n  of  most  features  of  interest. 

Provisions  are  available  to  link  HULL  to  the  Lagrangtan  TOODY  code  (REF 
2)  or  EPICS  <'ode  (REF  3)  for  calculations  in  two  and  ‘Free  dimensions  respec¬ 
tively. 

Development  of  HULL  (REF  4)  started  in  1971  at  the  Air  Force  Weapons 
Laboratory  (AfWL).  This  initial  version  was  bnii ..  to  solve  nuclear  weapons 
effects  problems  with  special  emphasis  on  nuclear  fireball  rise  ant  sir  bias* 
propagation.  It.  .1975  a  conventional  weapons  efiects  version  of  EUL;.  •  tEF  5) 
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was  developed  at  the  Air  Force  Armament  Laboratory  (AFATL).  This  later  ver¬ 
sion  oi  the  code  incorporated  new  physics  modules  for  calculating  hydrodynamic 
and  elastic/plastic  forces.  In  addition  the  code  architecture  was  modified  to 
permit  vector ization  on  pipe-line  machines  such  as  the  CRAY.  HULL  has  been 
continuously  developed  by  the  AFATL  and  AFtfL  since  1975  with  primary  manage¬ 
ment  by  AFATL.  Previous  versions  or  offshoots  of  HULL  are  not  currently 
supported  by  the  AFATL/ AFWL  team. 

The  extension  of  HULL  to  provide  c  link  to  the  Lagrangian  EPIC3  code  was 
sponsored  by  the  Army  Ballistic  Research  Laboratory  (BRL).  This  effort  was 
begun  in  September  1979.  It  is  the  purpose  of  this  report  to  document  this 
latest  version  of  HULL  (Version  113)  and  to  serve  as  a  user  guide  for  running 
HULL  as  an  element  in  linked  Eulerian/Lagrangian  penetration  calculations  in 
three  spatial  dimensions. 
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SECTION  II 
TEE  HULL  SYSTEM 


The  HULL  system  of  codes  is  composed  of: 


(1) 

(2) 

(3) 

(4) 

(5) 

(6) 


KEEL  -  an  initial  condition  (or  mesh)  generator, 

HULL  -  the  finite  difference  code, 

PULL  -  the  graphics  package, 

BOW  -  the  dnmp  tape  library, 

PLANK  -  the  alternate  input  generator  for  SAIL 

CONTROL  -  an  interactive  program  status 
monitor,  and 


(7)  SAIL  -  the  system  file  manager  and  executive 
proce  s  sor . 


All  programs  or  codes  in  the  HULL  system  including  SAIL  itself  are  main¬ 
tained  on  a  single  SAIL  program  library  file.  Compilation  of  any  program 
element  of  HULL  entails  execution  of  SAIL  to  process  the  program  library, 
select  the  proper  subroutines  and  options  desired  for  the  particular  applica¬ 
tion,  and  produ  -e  a  file  which  will  serve  as  input  to  the  compiler.  This 
processing  activity  is  done  to  keep  storage  requirements  and  central  processor 
time  requirements  to  the  practical  minimum  consistent  with  code  comprehension 
and  structure  flexibility.  The  use  of  SAIL  for  running  HULL  and  its  asso¬ 
ciated  service  programs  will  be  addressed  iD  this  document.  Other  documents 
which  describe  the  HULL  code  are  limited  to  specific  features  or  elements  of 
HULL  (REF  6,  7.  and  8).  Use  of  SAIL  for  program  alteration  and  file  main¬ 
tenance  is  covered  by  Reference  8. 


Figure  1  illustrates  the  interactions  of  the  program  elements  of  tne  HULL 
system.  Each  of  the  programs  are  described  more  completely  in  the  respective 
section  of  this  report. 


The  general  setup  and  running  of  a  HULL  calculation  (assuming  that  an 
executable  SAIL  module  and  the  HULL  system  file  are  available  on  the  computer 
system)  starts  with  the  formulation  of  an  input  deck  for  KEEL.  This  input 
deck  contains  all  the  specifications  for  formulating  the  initial  conditions 
for  the  calculation.  The  KEEL  input  is  first  read  by  PLANK.  PLANK  processes 
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the  deck  to  determine  the  mesh  size  and  options  to  be  nsed  for  this  run. 
FLANK  then  generates  a  disk  input  file  for  SAIL  which  details  the  structure  of 
program  KEEL.  SAIL  uses  this  disk  file  to  process  the  HULL  system  library 
file  and  assemble  the  SEEL  compiler  '.put.  The  compiler  input  is  compiled  and 
the  resulting  code  is  loaded  and  executed  as  program  KEEL.  KEEL  reads  the 
same  input  deck  used  for  the  generation  of  KEEL.  KEEL  then  produces  the  HULL 
problem  initial  conditions  and  puts  the  results  on  file  TAPE4.  TAPE4  is 
either  a  physical  tape  or  a  permanent  file  disk  device  which  will  store  the 
data  for  running  HULL  (or  PULL).  A  second  optional  tape  file  output  is 
produced  if  detailed  time  history  data  at  several  arbitrary  HULL  mesh  points 
are  desired.  This  is  called  the  station  tape  and  it  resides  on  logical  unit 
TAPE9.  The  station  tape  is  required  for  linking  HULL  to  the  Lagrangian  TOODY 
or  EPICS  codes. 

The  output  data  from  KEEL  on  TAPE4  is  used  as  input  to  PLANK  to  help  SAIL 
tailor  the  HULL  program  for  a  particular  run.  This  tailoring  process  can  be 
modified  by  the  HULL  input  dock  as  will  be  described  later.  PLANK  and  SAIL 
produce  the  compiler  input.  When  HULL  is  compiled  and  loaded,  it  reads  the 
initial  conditions  from  TAPE4  and  starts  the  calculation.  Snap-shots  of  the 
proceeding  calculation  at  pre-determined  times  (from  input)  are  put  on  TAPE4 
for  later  analysis  or  restart.  Station  da;a  is  also  accumulated  and  recorded 
on  TAPE9  if  it  is  present.  TAPE4  snap-shots  or  dumps  are  also  made  periodi¬ 
cally  every  60  central  processor  minutes  to  protect  against  unscheduled 
machine  failures.  Thus  TAPE4  and  TAPE9  store  a  continuous  record  of  the 
running  calculation  that  is  sufficient  to  restart  the  calculation  at  any  of 
the  available  TAPE4  dump  rimes. 

TAPE4  and  TAPE9  are  the  permanent  record  of  the  calculation  for  subse¬ 
quent  analysis  and  editing.  PLANK  uses  these  tape  as  and  a  PULL  input  deck 
to  instruct  SAIL  in  the  generation  of  the  PULL  c o-le  PULL  produces  an  ex¬ 
haustive  set  of  graphical  outputs  which  can  be  selected  £ ',r  problem  analysis. 
The  PULL  code,  after  being  generated  by  SAIL  and  compiled  bj  the  operating 
system,  reads  TAPE4  or  TAPE9  and  produces  graphics  output  for  the  selected 
plotting  devices. 
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SECTION  III 

BASIC  EQUATIONS  AND  SOLUTION  ALGORITHM 


The  HULL  code  solves  the  partial  differential  equations  of  continuum 
mechanics.  Explicit  terms  for  heat  conduction  and  viscous  effects  are  not 
included.  The  equations  solved  in  finite  difference  form  by  HULL  are: 


>  +  P»i,i  -  0 


p4j  -Tij,i 


pg; 


(1) 

(2) 


where : 


pB  -  (TyUjl.j  -  -pujgj 
Tlj  =  + 


)  =  material  density  (grn/cc) 

i  =  material  velocity  (cm/sec) 

« J 

'ij  =  ^ij  “  ^ijP»  stress  tensor, 

P  =  hydrostatic  pressure  (dynes/cm1) 

>2^.  =  stress  deviator  (dynes/cm1) 

=  gravitational  body  force  (cm/sec*) 

l  -  I  +  1/2  u.u.,  total  specific  energy, 
1  I 

I  =  internal  specific  energy  (ergs/gm) 


2Y2 

subject  to:  “  ,  Y  =  flow  stress. 


The  deviatoric  strain  rate  tensor  can  be  expressed  as: 


(3) 

(4) 


5ij  ■ "  <ui,j +  "j,i)/2 "  Vi.i/3- 


(5) 
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The  deviatoric  stress  tensor  is  linearly  related  to  the  strain  tensor  by: 


Sij  "  2Geij 

where  G  is  the  elastic  shear  modulus. 


(6)  • 


rJL  a  5  1 1  C  x  uy  iS  muucj 


the  Von  Mises  yield  criterion: 


JVi 

S,.  S..  <  — r-*  T  =  tensile  flow  stress. 


If  the  stress  deviators  are  not  contained  within  this  measure,  they  are 
brought  back  to  the  yield  surface  by: 


3Vij 


1/2 


(7) 


The  initiation  of  failure  can  be  determined  by  one  of  three  criteria: 

(1)  Maximum  principal  stress  can  be  used  as  a  crude  means  to  deter¬ 
mine  if  a  spall  threshold  or  failure  in  plane  strain  has  occurred.  The 
principal  stresses  are  defined  by  solving: 

11^-06^1=0.  (8) 

Failure  is  initiated  whenevor: 

Maximum  [o]^.  F*Y 

where  F  is  a  constant  between  2.5  and  4. 

(2)  Maximum  principal  strain  can  also  be  used  to  determine  the 
onset  of  failure  in  plane  stress  or  ductile  failure  by  solving  the  analog  of 
equation  (8)  for  the  maximum  principal  strain  as: 


I e i j  -  X6 ^ I  =  °  (9) 

and  assuming  failure  occurs  whenever 
Maximum  [X]  e^ 

where  e^  is  the  tensile  strain  to  failure  of  the  specific  material. 
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(3)  Triaxial  states  at  which  failure  can  occur  are  based  on  a  P/ Y 
model  (REF  6).  In  this  formulation  failure  is  initiated  when  the  maximum 
principal  strain  (X,  equation  9)  exceeds  a  failure  surface  X^  =  f(P/Y)  which 
is  determined  by  experiment.  Failure  in  HULL  is  initiated  by  inserting  a 
minimum  numerically  significant  void  in  the  failed  zone  and  allowing  the 
tensile  forces  to  relax.  Subsequent  tension  or  shear  of  this  region  is  not 
supported  (i.e.  the  zone  defines  a  slip  surface  and  tensile  stresses  are 
zeroed) . 

DIFFERENCE  method 

The  solution  algorithm  will  be  expressed  only  for  the  three  dimensional 
version  of  HULL.  For  two  dimensions  see  Reference  4. 

Equations  (1)  through  (3)  are  solved  on  a  finite  difference  mesh  composed 
of  discrete  spatial  intervals  x^,  y^,  z^  in  the  three  dimensional  Cartesian 
coordinate  system  having  coordinates  x,  y,  and  z.  The  solution  is  advanced 
explicitly  from  the  initial  conditions  by  discrete  time  intervals  Atn.  The 
solution  at  any  time  tn  for  any  variable  f(x,y,z,t)  in  the  solution  space  is 

defined  by  f.  ,  .  =  f(x.,y.,z.,t  )  where: 

n  l  j  ic  n 

i-l 

\  Ax. 

xi  =  Xo  +  1  tXm  +  ~1X 
m=i 

j-i 

y .  =  y  +  'S  Ay  + 
j  J  o  £.m  2 

m=i 

zk  =  zo  +  Ifm  + 
n 

.  n  ±iz 

t  =  t  +  )  Az„  +  — m 
o  L  «  m 
m=l 

Note  that  on  this  mesh  Ax^  =  x.^  Xj. 
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This  definition  holds  analogously  for  y,  z,  and  t. 


The  variables  that  are  defined  for  each  point  in  the  mesh  at  time  tn  a>nd 

location  x.,  y.,  z,  are: 

1  J  K 


P 


U 

V 

w 


I 

H 

M 


XX 


yy 

s 


xz 


yz 


=  pressure  (dynes/cm1) 

=  x  component  of  velocity  (cm/sec) 

=  y  component  of  velocity  (cm/sec) 

=  z  component  of  velocity  (cm/sec) 

=  internal  specific  energy  (ergs/gm) 

=  mass  (gms) 

=  mass  of  the  nth  species  (gms) 

=  volume  occup’ud  by  the  nth  species  (cm3) 

=  x  coordinate  deviatoric  stress  (dynes/cm*) 

=  y  coordinate  deviatoric  stress  (dynes/cm2) 

=  shear  stress  normal  to  z  coordinate  (dynes/cm2) 
=  shear  stress  normal  to  y  coordinate  (dynes/cm2) 
=  shear  stress  normal  to  x  coordinate  (dynes/cm2) 


Since  the  deviators  have  the  property 


S  =  -S 

ZZ  XX 


S  is  not  explicitly  stored  as  a  mesh  variable. 

ZZ 

If  we  define  the  central  difference  operator  oZ(f.  .  )  as: 

J/ 


5  (fi,j,k}  (fi+l/2,,i,k  ~  fi-l/2,i,k) 

Axi 

then  we  can  express  the  calculation  for  the  velocity  components  as: 


(10) 


un+l  =  un  +  ^ 
Pn 


6x(SD  ~  Pn+1/2)  +  5y(Sn  )  +  6z(s“  ) 
xx  xy  xz 


(11) 


=  vn  + 

Pn 


6y(s“  -  Pn+1/2)  +  6x(Sn  )  +  8z(Sn  ) 
yy  '  xy  yz 


(12) 
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n+1 


n  ,  At 

y  +  — 


where : 


6Z(-Sn  -  SD  -  Pn+1/2)  +  6s(Sn  )  +  6y(Sn  ) 
sx  yy  xz  yz 


n+1/ 2  _  n 

i+1/2  “  Pi+l/2  2  ^  i+1/2  "  Vl'i+1/2J 


_ _ _  __  At  J2-  cx ,  n  > 

/o  ~~  Pjxi  /o  o  /o  *  /o' 


and: 


pi+lPi  +  piPi+l 


i+1/2 


Pi  +  Pi+1 

pC;+i/2  =  min  (Pc;J>k  .pc;+1>j>k> 


n  _  Pitti  +  P  i+J  ai+l 

li+l/2  Pi  +  Pi4.i 


Solution  of  the  conservation  of  energy  equation  is 
„n+l 


e”  + 

n 

P 


5*  (_p  un+l/2  +  sn  +  s*  y«  +  gn  wn} 

xx  xy  xz 


+  6y  (~Pvn+1/2  +  Sn  vn+1/2  +  Sn  un  +  Sn  wn) 

yy  xy  yz 


,  az  /  »v  n+1  / 2  ,  cn  n+1/ 2  ,  cn  n  „n  n. 

+  o  (“Pw  +S  w  +S  u+S  v) 

zz  xz  yz 


where  the  boundary  value  velocities  at  time  n+1/2  are  given  by: 

n+1/ 2  _  n  At  ,r,Dn  v 

ui+l/2  "i+1/2  n+1/2  6  (Pi+l/2) 

pi+l/2 

with  similar  forms  for  v  and  w  with 


n  piui  +  pi+l  ui+l 


u. 


i+1/2  p.  +  pi+1 


(13) 
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n+1/2  n 


At  ex  n 


pi+l/2  “  pi+l/2  (1  "  2  6  ui+l/2) 


P i+1/2  =  Maximum  (Pi  »  Pi+i) 


Each  of  the  S  Q  terms  are  calculated  in  a  manner  similar  to  S  as: 

ap  xz 


S“1+i n  iitm 


p  .  ,  1  S  +  p  .  S 

in  XX.  1  xx.+1 

pi  +  pi+l 


^i+1/2  =  MIN(vfi'  vlW 


and  vf^  =  the  fraction  of  zone  i,  j,  k  which  is  occupied  by  a  solid  material. 
The  time  n+1  deviators  are  calculated  by: 


,n+l 

» 

=  sn 

+  2G(8x(un+1/2) 

-  D)At 

xx 

XX 

,n+l 

> 

=  sn 

t  2C(6y(vn+1/2) 

-  D)At 

yy 

yy 

,n+l 

» 

=  sn 

+ 

G(6X(vn)  +  6* 

(un))At 

xy 

xy 

,n+l 

i 

=  sn 

+ 

G(8x(wn)  +  8Z 

(u11)  )At 

xz 

xz 

,n+l 

» 

=  sn 

+ 

G(6y(wn)  +  8 2 

(vn))At 

yz 

yz 

where : 


rcx,  n+1/2.  cy,  n+1/2.  ,  £z,  n+1/2. , 

D  =  [6  (u  )  +  bJ (v  )  +  6  (w  )J/3. 

The  stresses  are  then  compared  with  the  yield  surface  definition  by: 


J  =  S2  +  S2  +  S  S  +  S2  +  S2  +  S2  <  Y3 / 3 
xx  yy  xx  yy  xy  xz  yz  -x 


If  this  criterion  is  not  met,  then  each  stress  component  is  reset  by: 

1/2 


S 


ap 


li 

3J 


$ 


ap 
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where  S  .  represents  each  of  the  elements  of  the  stress  tenser, 

cp 

Up  to  this  point  the  entire  calculation  has  been  in  a  Lngrangiaa  frame¬ 
work.  The  next  step  is  to  indicate  the  •volume  change  due  to  spherical  forces 
or  pressure  and  prepare  the  system  for  the  fluxing  calculation*  This  is  done 
by  solving  the  continuity  equation  by: 

Vn+1  =  Vn(.U2d  t5» 


where  V  is  the  cell  volume.  Tho  equation  of  state  is  then  called  and  the 
pressure  in  the  new  zone  is  calculated.  If  the  tone  contains  more  than  one 
material,  an  iteration  procedure  is  used  to  equilibrate  the  pressure  of  each 
material  species  in  the  zozie  to  a  com  non  value..  This  is  done  by  reparti- 
tioning  the  zone  among  the  materials  present. 


Ti  a  calculation  is  returned  to  a  fixed  (or  uniformly  translating)  grid  by 
fluxing  material  from  the  displaced  to  the  fixed  grid.  The  volume  to  be 
fluxed  is  defined  as  the  sura  of  the  transfer  volumes  in  each  of  the  three 
directions,  which  are  given  by: 


v  t  n+l  a  .  .  .  ..a+l.-n 

V i+1/2  "  xi+l/2  "  x i+1/2  Ayj  Azk  V  /V 

„  .  n+1  n  ,  .  ,  ,n+l /vn 

Vj+l/2  “  {yj+l/2  yj  +1/2  A“i  A\  V  /  J 

,,  ,  n+1  n  .  .  .  ,,c+l/t,n 

Vk+l/2  ~  7'k+l/2  k+l/2)  Ayj  V  /V 


where 


V. 


i+1/2 


“  volume  (in  cms)  to  be  transported 
between  cell  i»j,k  and  cell  i+l,j,k 


n+1 

i+1/2 


x , 


_  n  .  n+1  A. 

~  xi+l/2  T  Ci+l/2  At 


n  n 

a  _  *t  Llill 


i+1/2 


(14) 

(15) 

(16) 


i  9 


n+1  n  n+1  n 

n+1  _  tti  +  ui+l  ^i+1 

Ui+l/2  n  n 

P;  + 


with  analogous  definitions  for  j  and  k. 

The  transfer  volume  defined  by  eqs.  14-16  is  apportioned  among  the  avail¬ 
able  materials  in  the  donor  cell  in  accordance  with  the  diffusion  limiter 
algorithm.  This  algorithm  is  ad  hoc  but  has  been  formulated  over  the  years  by 
observation  of  various  numerical  experiments  and  has  been  found  to  be  modestly 
successful.  The  algorithm  uses  nearest  neighbor  information  in  an  attempt  to 
unmix  mixed  material  cells.  This  is  done  by  assigning  a  relative  transport 
weight  to  each  material  in  the  donor  cell,  and  then  normalizing  to  the  trans¬ 
fer  volume  given  in  equations  14-16.  Each  direction  (up  to  6  in  3D)  is 
handled  separably  subject  to  the  constraint  that  the  sum  of  the  exiting 
volumes  for  each  material  shall  not  exceed  the  initial  volume  of  that 
material . 


The  relative  transport  weight  for  each  material  n  is  defined  by 


VnR  -  VnU 
VnR  +  VnU 


+  1.0 


2V 


nR 


V  _  V  1 
nR  +  nU  j 


(17) 


where  V  _  =  fractional  volume  of  material  n  in  the  receiver  cell 
nR 

V  „  =  fractional  volume  of  material  n  in  the  upstream  cell. 
nU 

The  fractional  volume  is  defined  as  the  volume  of  material  n  divided  by  the 

total  cell  volume.  Two  exceptions  to  equation  17  are  noted.  If  there  is  no 

material  n  in  the  donor  cell  then  W  =0,  and  if  V  n  and  V  „  are  zero  then  W 

r  n  R  nil  n 


Finally  the  transport  fraction,  A  ,  for  each  material  n  in  a  particular 

n  m 

direction  m  is  obtained  by  normalizing  the  relative  transport  weights  by 
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(18) 


A 

nm 


'CIV 

Ann 


W 

n 


where  Vp  =  transfer  volume  given  by  equation  14,  15,  or  16 

V  =  volume  of  material  n  in  donor  cell, 
n 

The  transport  fractions  A  are  subject  to  the  constraint  tha;  the  sum  over 

nm 

all  directions  (six  in  3D)  shall  not  exceed  unity  for  each  material  n.  If  it 
does,  suitable  corrections  are  made. 

The  transport  fraction  applies  not  only  to  the  volume  but  also  to  the 
mass  and  internal  energy  for  material  n.  That  is,  the  fraction  Anm  of  each  of 
these  quantities  will  be  transferred  from  the  donor  to  the  receiver  cell.  The 
sum  over  n  for  each  direction  defines  the  total  flux  of  mass  and  internal 
energy. 

The  transport  methodology  will  be  illustrated  for  one  direction  (flux  to 
the  right).  The  remainine  directions  are  analogous.  Consider  cell  i,j,k  at 
the  end  of  the  Lagrangian  phase,  denoted  above  as  time  n+1.  For  each  material 
n  the  amount  of  mass  transporting  is 

^nR  ^n^  donor 

while  the  volume  transporting  is 

b/  =  (A  _  V  )  . 
n  nR  n  donor 

and  the  internal  energy  is 


6 (MI;  -  [A  „(MI>  ] . 

n  nR  n  donor 


The  total  mass  fluxing  (to  the  right)  is  therefore 


bM8  =  J  bM^ 
n 
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while  the  internal  energy  is 


6(MI)R  =  5  6<M  I  )L\ 
i*  n  n 

n 


The  total  energy  is 

6<ME)a  =  6(MI)K  4-  1/2  hl/P'  (ud*  +  vd51  +  wd*) 

where  d  indicates  donor  cell  velocities.  The  transporting  material  carries 
with  it  a  proportionate  amount  of  momentum,  thus 

6(Mu)R  =  u.&M8 
o 

6 ( Mv ) R  = 

5(Mw)R  =  v.hbP 

d 

With  most  options  HULL  transports  certain  histories  from  cell  to  cell.  These 
histories  essentially  transport  with  the  mass,  hence  for  the  kth  history  H 

6HR  =  H  , 

k  dk 

where  d  again  indicates  the  donor  cell. 

Similar  terms  are  obtained  for  the  other  directions  fore  (F)  and  above 
(A),  and  from  previous  calculations  for  left  (L),  aft  (T),  and  below  (B). 

Terms  such  as  the  mass  M.n+.^ ,  are  redefined  from  the  displaced 

1  9  J  /  & 

Lagrangian  mesh  to  the  original  fixed  (or  uniformly  translating)  Eulerian  mesh 
eM.nt*  as  follows.  For  matorial  n 

i»  J,k 

eMn+1  -  Mn+1  +  5^  +  6MT  +  6#  -  5M11  -  8#  -  5MA 
n  n  nnnnna 

eVn+1  =  ya+1  +  hV1  +  5VT  +  5VB  -  SV5  -  SVF  -  6VA 
n  n  nnnnnn 

°(MI)n+1  =  (MI)a+1  +  6(MI)L  +  8(MI)X  +  8(MI)B  -  8(MI)I{  -  6(MI)F  -  8(MI)A 
n  n  nnnnnn 


2  2 


The  velocity  is  defined  by  conserving  momentum  with 


eun+1  =  [un+1  Mn+1  +  8(Mu)L  +  8(Mu)T  +  S(Mu)B 

-  8(Mu)R  -  6(Mu)F  -  8(Mu)A]  /}  eM“+1 

n 

evG+1  =  [vn+l  Mn+1  +  6(Mv)L  +  s(Mv)T  +  6(Mv)B  _  6(Mv)R 

-8(Mv)F  -  8(Mv)A]  /  J  eMG+3 

n 

ewn+1  =  [wn+1  Mn+1  +  6(Mw)L  +  6 (Mw)T  +  8(Mw)B  -  8(Mw)R 

-8(Mw)F  -  6(Mw) A]  /  }  eMn+1 

L  n 

n 

The  new  internal  energy  is 

eIn+1  =  [In+1  Mn+1  +  8{ME)L  +  6 (ME)T  +  6(ME)B  -  6lME)R  -  6(ME)F  -  6 ( ME ) A ] 

,  V  e,,n+l  r/e  n+1.2  ,  ,e  n+l,1  ,  ,e  n+1,2, 

/  )  M  -  1/2  [(  u  )+(v  }+(w  )  J 

L  n 

n 

The  new  histories  are  defined  by 

CIl£+1  =  (Il£+1  Mn+1  +  6il^'  +  6Il£  +  SH,®  ~  6HR  -  SIl£  -  8Il£)/  ^  eM”+1 


It  should  be  noted  that  in  general 

(VttK,l'*1)  *  )  e(HI)n+1 

L  n 

n 

as  these  terms  a-e  defined  above.  This  is  due  to  a  dissipation  of  kinetic 
into  internal  energy.  Totul  energy  is  conserved  by  adjusting  the  individual 
energies. 
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SECTION  IV 
EQUATIONS  OF  STATE 

This  section  describes  equations  of  state  employed  within  SOS  “5  in  the 
HULL  code  in  versions  109  and  later.  The  equations  are  reviewed,  along  with 
the  assumed  material  properties.  The  equations  do  not  generally  attou  pt  to 
accurately  cover  all  regions  of  P,  V,  I  space.  They  were  developed  to  pi.ov'Hc 
valid  material  response  in  problems  of  interest  in  conventional  munitions 
research  -  shaped  charge  jet  formation,  self  forging  fragment  formation- 
metal-on- m  etal  penetrations,  bomb  penetrations  into  concrete,  and  other 
similar  problems.  This  emphasis  has  meant  for  example,  that  great  care  was 
given  to  modelling  solid  and  liquid  metal  behavior,  but  metal  vapor  behavior 
was  simplified  to  an  ideal  gas  formulation.  Problems  of  conventional  muni¬ 
tions  research  interest  do  not  usually  require  highly  accurate  metal  vapor 
response. 

The  HULL  equations  of  state  (EOS  6)  have  been  implemented  in  EP1C3  with 
some  minor  simplifications.  Thus,  this  report  can  serve  as  documentation  for 
many  of  the  EPIC3  equations  of  state  as  well  as  those  in  HULL. 

In  what  follows,  the  term  equations  of  state  will  refer  to  the  determina¬ 
tion  of  all  stress-strain  relations  including  deviatoric  as  well  as  hydrosta¬ 
tic  behavior. 

1.  Hefal  Equations  of  State 

The  metal  equation  of  state  relationships  are  based  on  Mie- Gruneisen  and 
gamma-lew  gas  pressure  formulations  and  the  von  Mises  yield  criterion.  The 
yield  strength  can  be  a  bi-linear  function  of  accumulated  plastic  strain  (work 
hardening)  and  a  tri-lmear  function  of  internal  energy  (thermal  softening). 
The  internal  energy  density  (energy  per  unit  mass)  at  onset  of  melt  and 
vaporization  are  increased  with  increasing  density.  Increasing  internal 
energy  is  not  allowed  to  iAcrease  pressure  during  the  melting  and  vaporization 
phase  changes.  We  have  made  no  provision  for  solid-solid  phase  changes  al¬ 
though  such  paase  changes  could  be  added  relatively  simply. 
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A.  The  Pressure  Equations 


The  basic  HULL  pressure  equation  for  a  solid  or  liquid  metaJ  is  the 
Mie-Gruneisen  equation: 

P  =  PH  (l-I>/2)  +  p(I-I0)  +  P0 

where : 


Py  is  the  Hugoniot  pressure  given  by: 
PH  -  Bl(i  +  BjP2  +  B3p»  if 


=  Blfl 


if  p<0 


and  p  is  the  excess  compression: 

P  =  P/P0  -  1 
p  is  current  density 

p  is  ambient  density 
o 

3P 

is  the  Gruneisen  parameter  (V  given  by: 


I  is  the  internal  energy  per  unit  mass 

I0  is  the  internal  energy  per  uiut  mr.ss  at  ambient 
density  and  pressure  (1  atmosphere) 

I>0  is  ambient  pressure  and  is  calculated  from 


Po 


l 1-013  x  1 06  d ynes/cm2 
I'o  "  ' 


Min  (I,  I0) 


This  allows  P0  to  be  1  atmosphere  when  I  =  I0  and  to  decrease  smoothly  to  zero 
as  J  drops  below  the  ambient  level  (!„)• 
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The  Hugoniot  coefficients  are  calculated  from: 


B*  “  poV 

B.  =  Bx  [1+2  (s-l)J 

Bj  -  B1  [2  (s-1)  +  3  (s-1)2] 


where  C  and  s  define  the  linear  shock  velocity  (U  )  vs  particle  velocity  (U  ) 
o  s  p 

relationship. 


The  assumption  of  a  linear  U  vs  U  relationship; 

s  p 

U  =  C  +  s  *  U 
so  p 

combined  with  the  Hugoniot  jump  conditions  actually  results  in  the  Hugoniot 
pressure  equation: 

PH  =  poCo  +  1)^1  ”  M(s-l)l 

Expanding  the  denominator  by  the  binomial  theorem,  simplifying  and  deleting 
all  terms  with  the  mul  .pliers  beyond  p5  yields  the  cocffients  Bx,  B2  and  B3. 
This  expansion  process  has  also  deleted  the  pole  in  the  P jj  formulation. 

The  Gruneisen  parameter  assumption  namely: 


F  =  rnpn/p 

o  o 

is  considered  reasonable  for  metals  (REF  1)  and  results  in  well  behaved  pres¬ 
sures  at  high  values  of  density.  The  other  often  used  assumption  that  does 
not  change  with  increasing  density  results  in  pressure  maxima  for  relatively 
small  density  increases. 


The  internal  energy  density,  I,  used  in  the  Mie-Gruneisen  pressure 
ion  is  not  alwa 
tionship  is  seen  below: 


calculation  is  not  always  the  zone’s  internal  energy  density,  I  .  The  rela 

z 
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where  is  the  internal  energy  density  at  the  onset  of  melt 

1^^^  is  the  energy  required  to  complete  fusion  or  the  melting  process 

1^  is  the  internal  energy  density  at  the  onset  of  vaporization 


IVAP  is 


the  energy  required  to 


complete  vaporization. 


As  the  diagram  indicates,  I  is  equal  to  1^  minus  internal  energy  lost  in  the 
fusion  and/or  vaporization  process.  Thus  pressure  due  solely  to  energy  addi¬ 
tion  is  not  allowed  to  increase  during  fusion  and  vaporization  until  these 
processes  are  complete.  Density  increases  will  still  provide  pressure  in¬ 
creases  through  the  first  term  in  the  equation. 


Metal  energy  at  the  onset  of  melt  increases  with  increasing  confine¬ 
ment.  Thus  the  1^  and  Iy  should  be  changed  as  pressure  increases.  HULL  uses 
the  GRAY  metal  (REF  6)  fit  for  1^  as  a  function  of  volumetric  strain.  The 
function  is: 


TM  =  IM0  (1  +  +  h*%) 

where : 
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a 


=  2  ro  -  2/3 
b  =  (ro  -  1/3)  (2  r0  +  1/3)  -  1 
n  =  l  -  pq/p 

Ijjq  =  internal  energy  density  at  onset  of  melt  at  ambient 
density 

We  increase  Iy  using  the  same  function,  i.e.. 


Ipyg  and  ly^  a^e  constant  regardless  of  the  zone's  compression. 

The  HULL  equation  of  state  transitions  to  a  gamma-law  gas  at  internal 
energies  above  sublimation  energy  and  whenever  the  pressure  so  computed  is 
higher  than  the  Mie-Gruneisen  pressure.  The  pressure  comparison  ensures  a 
smooth  transition  between  equations  without  actually  determining  and  identi¬ 
fying  intersections  of  the  two  regions  in  P,  V,  I  space.  The  gamma-law 
equation  is: 


P  =  (y  -  1)  P  (I  -  I$) 

where  I  is  the  internal  energy  density  for  sublimation  and  y  is  C^/C^.  The  I 
in  this  equation  is  the  zone's  internal  energy — i.e.,  it  is  not  reduced  by 

XFUS  “  ■£V4P* 

There  are  a  number  of  const' aints  applied  to  the  Mie-Gruneisen  pres¬ 
sure  to  more  physically  model  metal  behavior  in  extreme  energy  and  density 
states.  The  constraints  are  more  easily  explained  by  expressing  the  Mie- 
Gruneisen  equation  as: 

P  =  Fx  (p)  +  F2  (p)  *  (I  -  Io)  +  Po 

When  material  in  a  zone  expands  under  tension  it  should  not  be  allowed 
to  do  so  indefinitely.  Use  of  the  void  insertion  fracture  model  in  HULL  can 
prevent  extreme  expansion  states  from  occurring.  However,  this  model  requires 
data  which  is  often  unavailable.  In  these  cases  the  following  constraints 
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v/ ill  provide  some  physical  realism. 

The  product  ~  I  )  is  limited  if  p  is  less  than  -0.25.  The 

equation  employed  is: 

F,  (I  -  I  )  =  [F,  (I  -  I  ) ]  20  (p  +  0.4) /3 

where  the  term  in  brackets  is  the  product  calculated  without  constraints.  The 
equation  provides  a  linear  descent  to  a  value  of  zero  for  the  product  at  p  =  - 
0.4.  Because  of  an  IF  check,  the  equation  is  not  invoked  unless  p  is  less 
than  -0.25.  The  point  of  controlling  this  product  is  simply  that  metals  at 
such  high  expansions  cannot  support  compression.  The  zone  contains  metal 
fragments  and  such  fragments  cannot  support  compression  when  they  are  so 
highly  distended.  Without  this  control,  compressive  pressures  could  still 
exist  with  high  enough  I  values. 

The  Gruneisen  parameter  is  limited  so  as  not  to  exceed  2F  .  In  highly 
expanded  states  the  Gruneisen  parameter  grows  without  limit  as  density  de¬ 
creases.  Since  the  Gruneisen  parameter  is  not  meant  to  be  applied  in  expanded 
zones,  it  is  reasonable  to  control  this  growth.  This  control  is  somewhat 

redundant  with  tho  F„(I-I  )  control. 

Z  o 

If  the  calculated  Mie- Gruneisen  pressure  is  tensile  there  are  several 
controls  which  are  applied. 

If  the  pressure  is  less  than  P,IT.,  (1-0. Ip),  it  is  reset  to  P,IT.,  (1- 

MIN  MIN 

O.lp).  Normally  P  ^  2S  set  t0  s  very  large  negative  number  (-1  x  10  2  0)  so 
that  this  control  is  not  invoked.  However,  the  same  equation  of  state  rou¬ 
tines  are  used  for  media  such  as  water  which  can  support  no  tension.  In  such 
cases,  this  control  is  required.  The  factor  (1-0. Ip)  provides  a  variable 
pressure  as  p  changes.  The  product  P^^  (l-0.1p)  grows  slightly  more  negative 
as  p  decreases. 

If  the  Mie- Gruneise n  pressure  is  tensile  and  p  is  less  than  -0.2,  the 
pressure  is  reset  to: 

-20  (p  +  0.2)  (80  -  P)  +  p 


29 


where  P  is  the  pressure  after  the  P^^  check.  If  p  is  less  than  -0.25  the 
pressure  is  reset  to: 

100  (1  +  p) 

These  two  controls  work  together  to  limit  total  pressure  to  very  small  posi¬ 
tive  values  along  a  linear  ramp  from  P  =  0  at  p  =  —  1  to  P  =  100  dynes/cm  2  at  p 
=  0.  This  provides  essentially  zero  pressures  with  a  small  positive  slope  as 
density  increases. 

The  metal's  tensile  capability  is  limited  as  its  internal  energy 
density  approaches  melt  energy.  To  use  this  control  the  quantity: 

JMY  =  IM  +  0"5  ^‘US 

is  calculated.  This  is  also  the  quantity  used  to  limit  yield  strength  and 
will  be  discussed  later  in  this  section.  If  I  is  greater  than  0.6  Iyy»  the 
calculated  tensile  Mie- Gruneisen  pressure  is  multiplied  by  the  factor: 

F  =  (IMY-I  )(IMY_I)  if  IMY  >  °*61MY 
MI  O 

"  °  if  1  >  *MY 

and  then  reset  to 

F  ’  P  +  100  (1  +  p) 

This  control  limits  tensile  values  of 
phases  the  pressure  into  the  distended 
discussed. 

Compressive  Mie-Gruneisen  pressures  are  limited  to  the  region  of  20  x 
10  1  1  dynes/cm2  (20  megabars).  Areal  material  cannot  achieve  such  pressures 
without  also  having  internal  energies  considerably  greater  than  Is«  In  this 
case,  the  equation  of  state  will  have  switched  to  a  gamma-law  gas  which  has  no 
compressive  controls  on  it.  Therefore,  if  a  call  is  made  to  the  equation  '•f 
state  where  P  exceeds  20  megabars  it  is  reset  to: 


pressure  as  I  increases  to  Iuv  and 

M  X 

fragment-filled  zone  curve  previously 
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P  =  2  x  1013  +  2  x  10*3  (|i-Mlo)/(100  -  p20) 


where  p20  is  the  value  of  ji  when  P  is  2  x  10  13.  This  provides  a  gentle  slope 
with  density  for  values  of  p  greater  than  p20. 

Before  adjusting  yield  strength  the  equation  of  state  routines  provide 
a  few  more  variables  required  in  other  sections  of  HULL. 


The  mixed-zone  equation  of  state  iterator  (REF  5)  requires  £*— — )  , 

I 

9P 

where  x  is  the  specific  volume,  1/p.  HULL  calculates  (r — )  and  multiplies  by 

9p  j 

-pa  to  obtain  this  quantity. 

The  adiabatic  bulk  sound  speed  is  required  for  several  purposes, 

including  time  step  control.  In  the  absence  of  heat  conduction,  the  adiabatic 

sound  speed  is  the  isentropic  sound  speed,  C  ,  calculated  from  the  thermo- 

s 

dynamic  identity  (REF  1): 


C 

s 


(3£)  =  (IP) 

' 9p' S  V9p'l 


+  ^<jf>0 

p  91  p 


For  both  pressure  derivatives  HULL  takes  account  of  the  controls  invoked  and 
computes  derivatives  along  appropriate  paths. 

An  approximate  temperature  is  calculated,  for  edit  purposes  only,  from 
the  equation: 


TCK)  =  I/C  +  T 
v  o 

where  Cy  is  calculated  from: 

Cy  =  3R/AW 

7 

R  =  gas  constant  =  8.134  x  10 
AW  =  the  metal's  atomic  weight 

and  where  Tq  is  the  temperature  required  to  force  T  to  a  value  of  288°K  when  I 
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is  IQ. 


B.  Yield  Strength  Equations 


If  the  W  0  R  K  option  is  set  to  1  in  the  KEEL  (generator)  input  deck, 
yield  strength  is  varied  as  a  function  of  accumulated  plastic  strain  in  a  form 
as  seen  below. 


e 

p 


where  YLDST  is  the  initial  (zero  plastic  strain)  yield  strength  and  YLDMAX  is 
the  ultimate  yield  strength  which  occurs  when  =  EPLAST.  If  WORK  =  0,  the 
ultimate  yield  strength  is  used  for  all  values  of  a^. 


Tho  calculation  of  e  incremental  strains  and  stresses,  and  resetting 
of  doviatoric  stresses  exceeding  the  yield  limit  are  carried  out  in  the  usual 
fashion  and  completely  described  in  the  previous  section.  These  calculations 
will  not  bo  repeated  here. 


After  an  ambient  energy  (possibly  work  hardened)  yield  strength  is 
calculated,  it  is  multiplied  by  a  thermal  softening  factor  which  varies  from  1 
at  ambient  energy  to  0  at  I^y.  Ijjy  is  the  melt  energy  at  which  all  yield 
strength  disappears.  UULL  sets  this  value  to: 


IMY  IK  +  °*5IFUS 

where  1^  is  the  onset  of  melt  energy  increased  as  a  result  of  increased 
density  and  IriT0  is  the  energy  absorbed  during  the  melting  process.  This  is  a 
somewhat  arbitrary  assignment.  There  will  be  some  point  during  the  fusion 
process  when  enough  centers  of  melt  exist  to  drive  yield  strength  to  essen¬ 
tially  zero  values.  The  exact  position  of  this  point  is  unknown  and  will 
probably  never  be  known  since  observables  such  as  pressure  and  temperature  do 
not  vary.  It  does  not  seem  that  the  metal  would  have  to  be  entirely  melted  to 
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essentially  zero  any  strength  characteristics.  On  the  other  hand,  when  &  few 
very  small  centers  of  melt  are  being  created — i.e„  at  the  onset  of  melt — it 
does  not  seem  that  the  strength  of  the  metal  would  be  appreciably  affected. 
To  avoid  both  extr^-es  HULL  uses  the  midpoint  of  the  fusion  process. 


The  thermal  softening  factor  is  detemined  frcm  a  trilinear  curve  as 
seen  below  : 


THERMAL 

FACTOR 


where  the  pairs  (YFRACl.IFRACl) ,  ( YFRAC2 , IFRAC2 )  are  material  constants. 


HULL  employs  the  ratio  1/ Ijj  y  instead  of  temperature  to  calculate  a 
thermal  softening  factor  simply  because  the  relationship  between  temperature 
and  internal  energy  is  virtually  unknown  at  densities  other  than  ambient.  The 
HULL  technique,  in  effect,  implies  that  there  is  a  constant  relationship 
between  I/IMy  and  T/TJ(  at  any  given  density  but  does  not  attempt  to  determine 
the  value  of  the  constant. 


It  is  possible  that  a  highly  expanded  zone  could  have  a  significant 
value  of  yield  strength,  when,  in  reality,  the  rubble  in  the  zone  would  be  too 
distended  to  support  any  strength.  To  prevent  this,  yield  strength  is  multi¬ 
plied  by  another  factor  in  addition  to  work  hardening  and  thermal  softening. 
This  factor  is  1  for  p  values  greater  than  -0.2.  It  varies  linearly  from  1  at 
p  =  -0.2  to  0  at  p  =  -0.25  and  then  is  0  for  all  values  of  p  less  than  -0.25. 

HULL  does  not  calculate  an  explicit  multiplier  for  strain  rate  varia¬ 
tions  of  yield  strength.  Metals  demonstrate  the  property  of  having  an  almost 
constant  yield  strength  as  strain  rates  vary  from  10}  sec-1  to  10J  sec-1. 
Since  conventional  weapon  calculations  lie  almost  completely  within  this 
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range,  HULL  assumes  that  the  strength  values  used  will  be  appropriate  within 
this  range. 


HULL  has  built-in  values  for  all  metal  material  properties.  Strength 
properties  can  be  changed  through  use  of  the  c'ATPROP  option  in  the  HULL  input 
deck.  Other  property  changes  require  changing  DATA  statements  in  the  code. 
Tables  I  and  II  summarize  the  built-in  property  values  for  materials  employing 
the  Mie-Gruneisen  equation  of  state.  The  Hugoniot  data  was  obtained  from 
References  1,  10  and  11  and  is  considered  very  accurate.  Melt  and  vaporiza¬ 
tion  energy  values  were  obtained  from  many  sources  and  considered  reasonably 
accurate  unless  noted  otherwise  in  the  tables.  Yield  strength  values  and 
thermal  softening  data  were  also  obtained  from  many  sources  including  material 
property  tests  at  the  U.S.  Air  Force  Materials  Laboratory,  the  fitting  of 
calculational  results  to  ambient  and  elevated  temperature  cylinder  impact 
tests  and  the  fitting  of  calculational  results  to  self  forging  fragment  for¬ 
mation  experiments.  However,  users  should  consider  all  of  the  strength  data 
as  suspect  and  attempt  to  obtain  data  valid  for  the  particular  metal  alloy  of 
concern.  When  such  data  has  been  found,  or  generated,  the  user  can  override 
built-in  strength  properties  by  entering  the  following  cards  in  HULL  input: 


MATPROP 

MAT  =  N 

YLDST  =  Y1 

YLDMAX  =  Y2 

EPLAST  -  E 

ENDPROP 

YFRAC1  =  YFj 
YFRAC2  =  YF2 

IFRAC1  -  IFj 
IFRAC2  =  IF2 

where  Y^,  Y^,  E,  YF^,  IF^,  YF2<  an^  ~^2  are  new  values>  and  N  is  the 

material  number  in  the  HULL  calculation.  The  example  illustrates  all  stress 
surface  values  being  changed.  Not  all  values  have  to  be  changed,  and  more 
than  one  material  can  be  altered. 

2.  Plain  Concrete  and  Sand  Equations  of  State 

HULL  contains  a  plain  concrete  equation  of  state  model  developed  to  pro¬ 
vide  reasonable  response  in  the  1  to  300  kilobar  range.  It  is  based  on  a  very 
small  amount  of  data,  since  most  available  data  is  at  stress  levels  well  below 
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TABLE  I 


HULL  MATERIAL  PROPERTIES 


MATERIAL 

HULL 

NAME 

P 

(gm/cc) 

co 

(10scm/5  ec) 

s 

0 

PMIN 

(dynes /cm* 

Aluminum 

AL 

2.71 

5.38 

1.337 

2.1 

-l.xlO*0 

Ceramic 

<"203> 

CERAM 

3.9 

6.9 

1.45 

0.5 

- l.xlO*0 

Copper 

CU 

8.9 

3.958 

1.497 

2.0 

-l.xlO*0 

Iron 

FE 

7.86 

4.61 

1.73 

1.67 

-l.xlO*0 

Lead 

PB 

11.34 

2.092 

1.452 

2.0 

-l.xlO*0 

Nickel 

NI 

8.9 

4.8 

1.254 

1.83 

-1.X10*0 

Nylon 

NYLON 

1.14 

2.29 

1.63 

0.87 

-10x10* 

Plastic 

(polyrubber) 

PLAST 

1.012 

0.854 

1.865 

0.9 

-10x10* 

RIIA 

(thin  plate) 

RHA 

7.86 

4.61 

1.73 

1.67 

-l.xlO*0 

Stainless 

Steel 

SSTEEL 

SAME 

PROPERTIES 

AS  IRON 

Tantalum 

TA 

16.6 

3.423 

1.214 

1.69 

-l.xlO*0 

Teflon 

TEFLON 

2.16 

1.34 

1.93 

0.9 

-10x10* 

Tungsten 

W 

18.1 

4.0 

1.268 

1.58 

-l.xlO*0 

Uranium 

U 

19.05 

2.48 

1.53 

1 . 6 

-l.xlO*0 

U-0 .75%TI 

UTI 

18.9 

2.48 

1.53 

1.6 

-l.xlO*0 

Water 

WATER 

1.0 

1.483 

1.75 

0.28 

0.0 
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TABLE  I  (COKT) 


HULL  MATES I AL  PROPERTIES 


IM0 

IFUS 

JV0 

IVO+IVAP 

*s 

y-1 (vapor 

MATERIAL 

(10*ergs/gm) 

(10»ergs / gm) 

(109ergs/gra) 

(10’ergs/go) 

(10*ergs/gm) 

equation) 

Aluminum 

7.3 

3.96 

32.0 

100.0* 

119.2 

0.25 

Ceramic 

(Al2°3> 

9.55 

2.1 

50* 

100.0* 

100.0* 

0.25 

Copper 

4.78 

2.12 

16.3 

63.0 

53,1 

0.25 

Iron 

7.4 

2.74 

22.4 

86.8 

74.2 

0.25 

Lead 

0.67 

0.262 

2.96 

11.2 

9.4 

0.25 

Nickel 

6.6 

3.1 

20.1 

84.7 

72.9 

0.25 

Nylon 

1.9 

2.1 

50.0* 

100.0* 

100.0* 

0.25 

Plastic  20.0* 

(polyrubber) 

2.0* 

50.0* 

100.0* 

100.0* 

0.25 

RHA 

(thin  pla 

7.4 

te) 

2.74 

22.4 

86.8 

74.2 

0.25 

Stainless 

Steel 

SAME 

PROPERTIES  AS 

IRON 

Tantalum 

4.32 

1.59 

12.5 

48.3 

43.2 

0.25 

Teflon 

20,0* 

2.1* 

50.0* 

100.0* 

100.0* 

0.25 

Tungsten 

4.77 

1.84 

13.6 

58.4 

46.0 

0.2 

Uranium 

1.4 

0.49 

5.5 

28.2 

22.7 

0.25 

U-0.75%TI 

1.4 

0.49 

5.5 

28.2 

22.7 

0.25 

Water 

7.1 

3.35 

23.8 

44.7 

44.7 

0.25 

*NO  DATA,  ARBITRARY  VALUES  ASSIGNED. 
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TABLE  II 


HULL  STRENGTH  PROPERTIES 


SHEAR 


YLDST  WYLDST  MODULUS 

(10*  (10*  (10*1 


MATERIAL 

DINES/CM*) 

DYNES /CM*) 

EPLAST 

YF1 

IF1 

IF2 

EF2 

DYNES/CM* 

Aluminum 

2.9 

2.9 

0.0 

0.9 

0.5 

0.9 

0.5 

2.69 

Ceramic 
(AljO, ) 

80.0 

80.0 

0.0 

0.21 

0.45 

0.05 

0.75 

10.0 

Copper 

1.2 

4.5 

0.3 

0.78 

0.41 

0.78 

0.41 

4.64 

Iron 

4.69 

5.5 

0.3 

0.9 

0.5 

0.9 

0.5 

6.41 

Lead 

0.3 

0.3 

0.0 

0.9 

0.8 

0.9 

0.8 

1.113 

Nickel 

6.0 

6.0 

0.0 

0.8 

0.5 

0.8 

0.5 

8.9 

Nylon 

0.5 

0.5 

0.0 

0.9 

0.5 

0.9 

0.5 

0.368 

Plastic 

(Polyrubber) 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

RIIA 

15.0 

15.0 

0.0 

0.9 

0.5 

0.9 

0.5 

6.41 

Stainless 

Steel 

6.89 

10.0 

0.3 

0.9 

0.5 

0.9 

0.5 

7.3 

Tantalum 

5 

5 

0.0 

0.5 

0.4 

0.5 

0.4 

6.56 

Teflon 

0.5 

0.5 

0.0 

O.S 

0.5 

0.9 

0.5 

0.233 

Tungsten 

20.0 

20.0 

0.0 

0.9 

0.5 

0.9 

0.5 

16.0 

Uranium 

8 

15 

0.5 

0.9 

0.5 

0.9 

0.5 

5.1 

U-0 .75%TI 

8 

15 

0.5 

0.9 

0.5 

0.9 

0.5 

5.1 

Water 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

J  7 


1  kilobar.  However,  the  equation  of  state  ha®  been  used  successfully  in 
calculations  of  penetrations  into  concrete  at  velocities  from  a  few  hundred  to 
over  1,000  feet/second  (REF  12). 

The  model  was  constructed  primarily  from  Hugoniot  data  on  one  specific 
concrete  (REF  13)  and  static  yield  strength  data  on  s,veral  concretes  (REF 
14).  Assumptions  were  made  for  a  value  for  Poisson's  ratio,  a  shear  modulus 
formula,  a  plastic  flow  rule  and  other  parameters. 

Gregson  (REF  13)  generated  Hugoniot  data  to  stress  levels  over  500  kilo- 
bars  for  a  granite  aggregate  concrete.  The  aggregate  had  an  average  diameter 
of  1/8  inch.  The  initial  density  varied  from  mix  to  mix  but  averaged  2.185 
gffi/cc.  The  solid  density  had  an  average  of  2.673  gm/cr.  The  concrete  con¬ 
sisted  of  18  percent  voids,  25  percent  granite  aggregate,  20  to  25  percent 
quartz  grains  and  25  to  30  percent  cement  paste  by  volume.  The  small  aggre¬ 
gate  size  was  dictated  by  his  shock-generating  method,  which  consisted  of  gas- 
gun  launched  flyer  plates.  He  employed  several  measurement  techniques  to 
establish  stress  and  velocity.  The  final  data  converted  to  stress  vs  excess 
compression,  p,  appears  in  Figures  2  and  3. 

The  lowest  point  in  Gregson's  data  is  at  0.75  kilobars.  This  is  a  fairly 
consistent  precursor  value  seen  In  his  stress  vs  time  data.  HULL  uses  that 
point  to  define  the  density  ratio  at  which  significant  cracking  and  void 
filling  occurs. 

The  estimated  unloading  path  on  the  figures  is  simply  one  which  allows  the 
concrete  to  return  to  its  solid  density  when  loadod  past  the  point  where  all 
voids  are  filled.  It  is  a  straight  line  which  connects  with  the  loading  curve 
at  the  total  void  filled  point  with  u  lope  equal  to  the  loading  curve  at  tha 
point.  There  is  no  data  for  this  line  beyond  the  solid  density  point  at  zero 
stress. 

The  yield  strength  of  concrete  increases  as  confining  pressure  is  in¬ 
creased.  The  only  data  available  at  high  confining  pressures  is  static  data 
generated  by  Chinn  (REF  14).  Chinn's  results  are  shown  in  Figures  4  and  8 
along  with  some  typical  low  stress  level  results  from  McHenry  (REF  15).  The 
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yield  strength  and  pressure  are  normalized  by  f',  the  concrete's  unconfined 

c 

compressive  strength.  Two  theoretical  points  are  shown  on  Figure  4.  One 

defines  failure  (flow)  at  the  unconfined  compressive  strength  point,  i.e.,  Y  = 

f  at  P  =  1/3  f'.  The  other  is  included  because  a  minimum  pressure  criterion 
c  c 

is  used  to  define  failure.  This  minimum  pressure  was  chosen  as  that  occurring 
in  an  unconfined  tensile  test  assuming  failure  when  stress  is  0.1  f'  (REF  16). 
It  is  numerically  pleasing  to  set  Y  =  0  at  this  failure  point. 

The  yield  strength  curves  are  approximated  in  HULL  with  three  straight 
line  segments: 

Y  =  0  if  P  i  -O.lf' 

c 

Y  =  3  (P  +  O.lf ' /3) /I .1  if  -0.1f'/3  i  P  i  f'/3 

c  c  c 

Y  =  P  +  2/3  f'  if  f'/3  <  P  i.  30f ' 

c  c  c 

Y  =  30.67  f'  if  P  >  30  f' 

c  c 

the  limit  at  P  =  30f'  is  not  indicated  in  Chinn's  data.  It  is  included 

c 

because  it  is  felt  there  should  be  a  limiting  value.  It  is  included  at  the 
level  shown  simply  because  it  is  not  known  to  exist  at  a  lower  level. 

It  should  be  emphasized  that  the  yield  strength  data  was  developed  in 
static  tests.  HifLer  yield  strengths  may  exist  at  strain  rates  of  interest  in 
conventional  munitions  work,  but  data  simply  does  not  exist. 

Using  the  previous  equations  for  Y,  Gregson's  Hugoniot  stress  data  was 
reduced  to  provide  an  Hugoniot  pressure  vs  excess  compression  curve  by  simply 
subtracting  2/3  Y  from  the  stress  data  assuming  an  f^  value  of  5,000  PSI. 
Figure  6  is  a  plot  of  this  curve  to  166  kilobars. 

Below  the  precursor  pressure  level  of  0.358  kilobars  (1/3  x  0.75  kilo- 
bars),  Hugoniot  pressure  is  fit  with  a  straight  line: 

PH  =  V 

where  Kq  is  determined  from  the  empirical  Young's  modulus  relationship  (REF 

10) 
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E  =  33W1,S  (r )0,5PSI  (W  in  lb/ft3  =  density) 
and  the  relationship 

K  =  E/3(l-2cr) 

with  a  value  of  Poisson's  ratio,  a,  chosen  as  0.2. 


Designating  the  value  of  p  at  0.358  kiiobars  as  p^, 
pressure  was  fit  with  the  following  equations: 


the  loading  Hugoniot 


PH(Kb) 

V 

0.358  +  78.62  (p-p  ) 

8  +  130  (p  -  0.1) 

21  +  420  (p  -  0.2) 

63  +  650  (p  -  0.3) 

76  +  625.5  (p  -  0.32)  % 

+  1720.19  (p  -  0.32) 

150  +  39.68  (p  -  0.414; 

166  +  360.5  (p  -  0.572) 

+  864  (p  -  0.572) 

The  maximum  density  unloading  curv 


RANGE  OF  p 

-0  to  pE 
PE  to  0.1 
0.1  to  0.2 
0.2  to  0.3 

0.3  to  0.366  (lockup  point) 
0.366  to  0.414 

0.414  to  0.54  (solid 

-  solid  phase  transition 

0.54  to  ® 

was  fit  by  the  equation: 


P„  =  784  (p  -  0.223) 

II 

HULL  maintains  a  zone  variable,  p^,  for  concrete  which  defines  the  maximum 
value  of  p  previously  seen  by  the  zone.  For  values  of  p^  less  than  p^, 
loading,  unloading  and  reloading  all  occur  along  the  same  elastic  path.  For 
values  of  pM  greater  than  p  at  lockup  (0.366)  loading,  unloading  and  reloading 
take  place  along  the  same  path.  For  p^  valuos  between  these  extremes,  un¬ 
loading  and  reloading  up  to  the  virgin  loading  curve  take  place  along  K  lines 

r 

where  K  is  chosen  to  vary  linearly  between  an(j  784  kiiobars. 
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The  shear  modulus  is  calculated  assuming  that  Hooke's  law  holds.  That  i*s: 


3(l-2or) 

G  =  - 

2(l+cr) 


Deviators  are  reset  as  though  the  yield  surface  were  a  von  Mises  surface  — 
i.e.,  not  a  function  of  pressure.  The  resulting  plastic  flow  rule  is  incon¬ 
sistent  with  the  Mohr-Coulomb  yield  surface  proposed  for  the  concrete.  How¬ 
ever,  the  alternative,  such  as  a  moving  CAP  model  (REF  2),  is  computationally 
complex  and  cannot  be  justified  considering  the  source  of  the  yield  strength 
model.  Strain  rate  effects  are  probably  far  more  important  than  a  consistent 
flow  rule. 

The  Hugoniot  pressure  previously  discussed  could  be  combined  with  a 
Gruneisen  parameter  into  a  complete  equation  of  state.  However,  Gruneisen 
parameter  measurements  for  concrete  are  typically  around  0.1.  For  this  rea¬ 
son,  the  equation  of  state  ignores  the  Gruneisen  parameter  and  Hugoniot  pres¬ 
sure  is  used  as  total  pressure.  This,  of  course,  means  that  internal  energy 
changes  have  no  effect  on  pressure  and  the  equation  of  state  could  not  be  used 
for  energy  deposition  problems. 

The  variable,  f',  is  set  to  the  CGS  equivalent  of  5,000  PSI.  If  other 

strengths  are  desired,  it  can  be  changed  in  the  data  statement  that  now  sets 

it.  Changing  f'  will  vary  the  fracture  point  and  the  yield  surface.  It  will 
c 

have  no  effect  on  the  pressure. 

The  sand  equation  of  state  was  developed  to  simulate  a  dry  40  percent 
porosity  quartz  sand.  It  combines  low  pressure  hydrostatic  data  from  Refer¬ 
ence  14  with  pure  quartz  Hugoniot  data  from  Van  Thiel  (REF  10).  The  sand 
model  assumes  that  strength  is  negligible  at  pressure  levels  of  interest. 
This  is  a  valid  assumption  as  long  as  strength  levels  of  a  few  hundred  PSI  are 
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unimportant  in  a  calculation.  There  is  evidence  that  sand  can  achieve  such 
small  levels  of  strength  when  highly  crushed. 

Sand  has  a  very  small  Gruneisen  parameter  which  can  be  ignored.  The 
equation  of  state,  then,  is  simply  a  pressure  vs  excess  compression  fit  as 
seen  in  Figure  7. 


The  virgi; 


1  rs  „  -  -  - 

j.v/auxu.5  V  UX  V  V 


jy  «  y  n  **y  c  m  oil 


<5*u  vnt  1 


P  =  K  p 
o 

from  p  =  0  to  p  =  0.0084.  Kq  is  computed  from  the  initial  density  of  1.6 
gm/cc  and  the  initial  sound  speed  of  0.61  x  10Jcm/sec.  The  pressure  level  at 
the  end  of  the  elastic  point  is  then  0.05  kilobars  (approximately  50  atmos¬ 
pheres)  the  elastic  segment  is  followed  by  two  linear  crushing  regions: 

P  =  0.05  +  4.96p  (Kilobars) 

for  p  from  0.0084  to  0.2 
and  P  =  1.0  +  73. 3p  (Kilobars) 

for  p  from  0.2  to  0.425 

At  a  p  of  0.425  and  a  pressure  of  17.5  kilobars  the  sand  is  assumed  to  have 
crushed  out  all  voids  and  the  loading  curve  follows  the  quartz  Hugoniot: 

P  =  2 . 9 8ri / (l-2.33q)a  (Kilobars) 

where : 

q  =  1  -  pQ/p 

To  avoid  the  pole  in  this  Hugoniot,  the  Hugoniot  continues  as  a  constant 
modulus  curve  at  q  =  0.98/s. 

Unloading  and  reloading  below  the  lockup  point  of  17.5  kilobars  and  above 
the  elastic  point  of  0.05  kilobars  follows  a  straight  line  curve  which  has  a 
bulk  modulus  continuously  varying  from  that  at  the  lockup  point  to  that  along 
the  elastic  curve. 
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The  sand  is  not  allowed  to  sustain  any  tension. 


Both  concrete  and  sand  require  an  extra  zone  variable,  the  maximum  value 
of  p  ever  seen  by  the  material.  To  successfully  use  these  models  the  user 
must  set  NHIST  =  1  in  tne  KEEL  input  deck.  This  will  cause  PLANK  to  increment 
NHIST  by  one  which  will  allot  one  extra  variable  per  zone  after  the  variables 
required  for  carrying  stresses.  This  will  hold  ^or  concrete  and/or  sand. 


3.  Miscellaneous  Solid/Liquid  Equations  of  State 

There  are  a  few  equations  of  state  in  HULL  other  than  those  following  the 
Mie-Gruneisen  formulation  (the  materials  in  Tables  I  and  II)  and  the  concrete 
and  sand  models.  These  will  now  be  briefly  discussed. 

Undetonated  AFX-108  is  very  simply  modelled  for  use  in  low  velocity  explo¬ 
sive-filled  bomb  and  warhead  impacts.  The  material  name  in  HULL  is  AFX1.  The 
Hugoniot  pressure  is  calculated  from: 


F 


where  is  set  to  10.3  x  10*  dynes/cm  1  .  The  Gruneisen  parameter  is  set  to  a 
constant  value  of  1.1  and  pressure  is  computed  from  the  Mie-Gruneisen  equation 
with  no  controls  or  tensile  values.  The  only  control  is  a  compressive  one.  p 
is  not  allowed  to  exceed  0.?.  This  strange  equation  of  state  is  the  result  of 
having  only  a  single  loading  modulus  for  the  material. 


Tillotson  equations  of  state  exist  in  HULL  for  granite  and  tuff.  They 
were  developed  for  nuclear  cratering  calculations  and  should  not  be  considered 
valid  at  low  conventional  weapons  research  pressure  levels.  They  have  not 
been  exercised  with  FLUXER  =  3  and  the  user  is  advised  to  avoid  them  if 
possible. 


4.  Air  Equations  of  State 

Two  air  equations  of  state  are  implemented  in  EOS  6.  If  the  user  simply 
asks  for  AIR,  the  equation  of  state  uses  the  strange  fovm  : 


where  y  =  1.4  and  Iq  is  set  to  return  one  atmosphere  pressure  at  p  =  pQ.  The 
purpose  of  this  form  is  to  provide  air  at  approximately  the  correct  density 
but  with  a  very  low  sound  speed.  Small  amounts  of  air  are  occasionally 
trapped  in  zones  as  the  result  of  motion  of  more  substantial  media.  High 
pressures  in  this  trapped  air  can  drive  the  time  step  to  very  small  values  and 
dominate  the  problem,  when  the  pressure  in  the  air,  and  perhaps  even  the 
presence  of  the  air  itself,  is  unimportant. 

If  correct  air  pressures  are  important  then  the  user  must  set  AIRE0S=1  as 
a  SAIL  option.  AIRE0S=1  insures  that  the  DO  AN- NICK  EL  variable  gamma  air  model 
will  be  used.  This  is  a  highly  accurate  model  over  broad  ranges  of  density 
and  energy  and  is  described  in  detail  in  Reference  26. 

Because  the  code  was  used  to  simulate  nuclear  fireball  calculations  at 
various  heights  above  sea  level,  HULL  is  capable  of  simulating  several  atmos¬ 
pheres  (air  initial  conditions).  These  are  described  in  Reference  27.  If 
ATMOS  is  not  specified,  a  constant,  sea  level  atmosphere  will  be  assumed. 

5.  High  Explosive  Equation  of  State 

There  are  several  high  explosive  equations  of  state  used  in  HULL.  The 
detonation  products  of  standard  military  explosives  are,  for  the  most  part, 
represented  by  the  JWL  formulation  (REF  18).  In  this  formulation,  pressure  is 
computed  from: 


P  =  A  (1  -  "a)  e  "(  1/q)  +  B  (1  -  ~|)  e  (  2/q)  +  wpl 

R1  R2 

where  q  =  p/p  ,  I  is  internal  energy  per  unit  mass  and  A,  B,  Rx,  R  a  and  w  are 
material  constants. 

This  equation  of  state  is  implemented  for  Composition  2,  Octol  (REF  22) 
and  C4  (REF  23).  JWL  data  exists  for  implementation  for  several  more  explo¬ 
sives  (REF  18-23). 

The  JWL  equation  has  proven  especially  useful  in  accurately  predicting 
metal  acceleration.  Its  utility  largely  derives  from  the  manner  in  >  hich  the 
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material  constants  are  determined.  They  are  forced  to  fit  several  consistency 
criteria  suc’i  as  total  deliverable  energy  equal  to  that  measured  calorimetri- 
cally,  pressure  at  Chapm  an-Jouge t  (CJ)  density  equal  to  measured  CJ  pressures, 

etc.  The  remaining  parameters  (Rlf  Ra,  and  o>)  are  fit  through  an  iterative 
experim ental/calcuiational  procedure  involving  measured  velocities  of  cylin¬ 
ders  and  spheres. 

At  low  densities  .be  last  term  in  the  JWL  equation  dominates.  HULL  takes 
advantage  of  the  calculational  simplification  offered  by  this  single  term  at 
high  expansions.  Whenever  density  of  the  detonation  products  decreases  to 
one-tenth  initial  density,  pressure  is  computed  from  the  single  pwl  term. 

For  densities  greater  than  one-tenth  ambient  the  JWL  pressure  does  not 
proceed  to  zero  as  energy  proceeds  to  zero.  This  situation  is  avoided  in  HULL 
by  decreasing  pressure  whenever  energy  is  less  than  1  x  10*  ergs/gm.  The 
calculated  pressure  is  multiplied  by  1/1  x  10*.  This  arbitrary  factor  allows 
low  energy,  relatively  high  density  gases  to  react  more  realistically.  By  the 
time  energy  has  dropped  from  initial  4  to  5  x  10  10  erg/gm  levels  to  1  x  10*, 
the  explosive  has  completely  delivered  its  energy  to  surrounding  media  and  is 
merely  expanding  to  fill  in  the  available  grid  regions. 

The  JWL  equation  of  state  calculates  a  temperature  for  editing  purposes 
only  from  the  equation: 

T(°K)  =  C^I  +  C2 
~8 

where  is  1.65  x  10  and  is  1375.  These  constants  were  found  to  fit 
available  data  for  Composition  B  and  are  used  for  all  JWL  explosives. 

The  JWL  constants  A,  B,  R^,  and  w  are  designated  as  B^,  B^,  B^,  B^  and 
Bj  respectively  in  HULL,  Values  of  these  constants  are  listed  in  Table  III. 

HULL  can  provide  burn  progression  (in  2D  only)  based  on  the  CJ  condition 
(D  =  C+U)  and  the  direction  of  pressure  gradients  in  the  unburned  explosive 
(BURN=1)  or  a  programmed  burn  can  be  used  (BURN=2).  A  programmed  burn  is 
suggested  when  an  accurate  detonation  velocity  is  desired  and  when  zone  sizes 
are  such  that  pressure  at  the  burn  front  is  significantly  less  than  C^  pres¬ 
sure.  Values  of  detonation  velocity,  D,  are  included  in  Table  III  for  use  in 
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Bt!RN=2.  Pressure,  and  Y  qj  =  C^/Cv  values  at  the  CJ  front  are  a." so  included  in 
the  table  for  calculational  comparison  purposes. 

If  BURN1^  is  desired,  the  user  inserts  values  for  TDET  (time  that  the 
explosive  region  begins  burning),  VDET  (detonation  front  velocity),  and  XDET, 
YDET  and  ZDET  (ignition  point).  These  values  are  inserted  as  the  last  element 
in  the  PACKAGE  setting  up  the  explosive  region. 


TABLE  III 

J1L  EXPLOSIVE  CONSTANTS  (OGS  UNITS) 


EXPLOSIVE 

HULL 

NAME 

po 

ro 

A 

B 

R1  R2  w 

COMP  B 

(64  RDX,  36  TNT) 

CBBRN 

1.72 

4.92  x 
1010 

5.242  x 
101J 

0.07678  x 
10ia 

4.2  1.1  0.34 

OCTOL 

OCTBRN 

1.82 

5.27  x 
1010 

7.486  x 
1011 

0.1338  x 

10  ** 

4.5  1.2  0.375 

C4  C4BRN 

(91  RDX,  5.3 
DI(Z-ETHYLHEXYL), 

2.1  POLYISOBUTYLENE, 

1.6  MOTOR  OIL). 

1.601 

5.62  x 
1010 

6.0977  x 
10ia 

0.1295  x 
1011 

4.5  1.4  0.25 

EXPLOSIVE 

NAME 

D 

C 

Ycj 

P  . 

Cl 

COMP  B 

CBBRN 

7.98 

x  10* 

0.01082  x 

1011  2.706 

295  x  10» 

OCTOL 

OCTBRN 

8.48 

X  10* 

0.01167  X 

1011  2.83 

342  X  10» 

C4 

C4  BRN 

8.193  x  10* 

0.01043  x 

10la  2.838 

280  x  10» 

It  is  sometimes  convenient  to 

in  expanding  from  the  CJ  point  to 

calculate  the  work  performed  by 

som  e  density  p  . 

an  explosive 

The  JWL  adiabat  (isentr-ope)  is  given  by: 
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p  = 


A  e  "(RlV  + 


a  ~(Vn>  + 


C(n)“+1 


the  values  of  C  for  this  adiabat  are  also  given  in  the  table.  Integration  of 
dp/p2  from  prT  to  the  density  of  interest  will  provide  deliverable  energy. 

Other  detonation  product  equations  of  state  exist  in  HULL.  They  can  be 
used  (cautiously)  for  ANFO,  BARATOL,  TNT,  PBX9404  and  Pentolite.  They  use 
various  forms  for  the  equation  of  state  and  were  developed  for  nuclear  explo¬ 
sive  simulation  purposes. 


ANFO,  TNT,  and  Pentolite  detonation  products  pressures  were  fit  to  the 
LSZK  equation  (REF  24  and  25) 

P  =  BlPI  +  B2pB3 

with  the  parameter  values  as  seen  in  Table  III.  The  TNT  and  Pentolite  data 
are  considered  to  be  very  accurate  at  large  expansions  and  were  extensively 
tested.  The  equations  have  not  been  checked  for  near  C^.  behavior  accuracy 
which  would  be  dominant  in  metal  acceleration  problems.  We  suggest  applica¬ 
tion  of  the  JWL  equation  for  this  region.  The  ANFO  data  used  in  HULL  was  from 
a  preliminary  ANFO  model  developed  by  Lawrence  Livermore  Laboratory.  Better 
ANFO  data  now  exists  and  is  contained  in  Referunce  20.  HULL  has  not  been 
updated  because  of  the  lack  of  interest  in  ANFO  in  conventional  munitions 
work. 

The  Baratol  equation  of  state  uses  a  gamma  law  formulation  simply  because 
better  data  doe*  not  seem  to  be  available.  The  material  parameters  for  this 
explosive  are  also  given  in  Table  III.  Based  on  comparisons  with  other  explo¬ 
sives  we  can  state  that  the  constant  gemma  assumption  will  provide  metal 
velocities  close  to  20  percent  too  high. 

A  PBX9404  fit  to  the  JWL  equation  of  state  is  included  in  HULL.  The  fit 
was  further  modified  to  reproduce  airblast  data  in  the  range  of  15  kilobars 
and  below.  It  should  not  be  used  for  metal  acceleration. 

The  predetonation  equations  of  state  in  HULL  mostly  use  the  form: 
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P  =  +  B^p2  +  BgH3  +  apl 


This  provides  a  simple  approximation  to  pressure  for  solids  expected  to  load 
along  a  Hugoniot  and  for  which  subsequent  unloading  and  reloading  are  unim¬ 
portant.  The  pressure  is  not  allowed  to  become  negative  in  this  equation  of 
state. 

Temperature  is  calculated  for  edit  purposes  only  and  is  fit  using  the 
equation : 

T  =  C1I  +  C2  if  T>288°K 
T  =  Cgl  if  T<2S8°K 

Values  of  the  coefficients  used  in  the  pressure  and  temperature  calcula¬ 
tions  are  included  in  Table  IV. 


TABLE  TV 

PRKDKTON  AT I ON  EXPLOSIVE  PARAMETERS  (CGS  UNITS) 

HULL 


EXPLOSIVE 

NAME 

B1 

*2 

„bjl_ 

C1 

c2 

C? 

COM?  B 

COMP  B 

1.282 

X 

0 

1.1932  x 

3.36  x 

200 

1.101 

10 

X 

1011 

10*2 

10"8 

OCTOL 

OCTOL 

1.282 

1021 

X 

0 

1.1932  x 
1012 

3.36  x 

10  8 

200 

1.101 

10 

X 

TOT 

TOT 

1.282 

X 

0 

1.1932  x 

3.312  x 

200 

1.085 
10  ' 

X 

1011 

1012 

10 

PENTOLITE 

PEN 

1.319 

X 

3.614  x 

2.295  x 

3.312  x 
10~8 

200 

1.085 

10 

X 

1011 

1011 

1011 

ANFO 

ANFO 

2.4 

X 

1.08  x 

7.2  x 

8  *-RX 
10  8 

200 

1.081 

10 

X 

10A0 

1011 

1010 
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The  Baratol  (HULL  name  BARTOL)  predetonation  equation  is  a  Mie-Gruneisen 
form.  The  coefficients  are: 

C  =  2.4  x  10s 
o 

s  =1.66 

F  =0.2 

o 

Temperature  fits  the  form: 

T  -  V  +  c2 

where  =  3.36  x  10  ^  and  C2  =  200. 

6.  High  Pressure  Geologic  Equation  of  State 

A  high  pressure  equation  of  state  for  geologic  materials  was  added  to  HULL 
to  provide  the  necessary  material  properties  for  the  calculation  of  meteor 
craters.  It  is  designed  primarily  to  provide  a  good  fit  to  Hugoniot  data 
between  50  kilobars  and  1  megabar.  Many  details  of  the  Hugoniot  are  only 
grossly  approximated  below  50  kb  and  so  the  equation  of  state  should  be  used 
with  caution  in  this  area.  It  is  based  on  the  information  in  References  28 
and  29,  modified  to  allow  a  fit  to  a  Mie-Gruneisen  equation  of  state. 

The  Hugoniot  is  modelled  in  several  connected  sections.  The  Hugoniot 
lock-up  curve  is  the  curve  which  describes  material  behavior  for  a  geologic 
material  with  no  voids.  The  crush-up  curve  describes  how  a  material  with  a 
non-zero  void  fraction  behaves  prior  to  reaching  a  zero-void  condition.  The 
two  Hugoniot  curves  are  seen  in  Figure  8.  In  this  figure,  is  the  density 
at  ambient  conditions  and  p^.  is  the  density  at  ambient  pressure  for  the  zero- 
void  material.  The  crush-up  loading  curve  is  described  by  a  single  bulk 
modulus,  Kj.  The  lock-up  loading  curve  is  described  in  several  segments.  A 
constant  bulk  modulus  segment  connects  the  1-atmosphere  point  vith  the  point 
at  which  the  crush-up  curve  meets  the  lock-up  curve  (jij ,  P^).  This  constant 
bulk  modulus  curve  then  turns  into  a  variable  modulus  curve  between  and 
Pjjj.  The  curve  is  described  by  the  Schuster-Ise nberg  form 

PI1  -  Vr  -  <W  1“ 
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where  p  =  p/p  -  1  and  K  and  K  are  initial  and  maximum  bulk  moduli  and  u* 
r  r  o  m  r 

is  a  parameter  which  determines  the  density  range  between  the  points  at  which 
Kq  and  Km  dominate.  At  the  point  (pgj,  Pgj)  the  Hugoniot  shifts  again  to  a 
constant  modulus  form 


PH  “  KHI  +  PEI 

the  unloading  path  taken  depends  upon  p  ,  which  is  the  maximum  value  of  u 

bi  ax  r 

seen  by  the  material.  If  Pmax  is  greater  than  then  all  unloading  and 

reloading  occurs  along  the  single  lock-up  curve.  If  u  is  less  than  u.  , 

then  unloading  and  reloading  occur  along  a  path  parallel  to  the  lock-up  curve. 

The  Hugoniot  pressure  is  unlimited  in  tensile  extent  but  the  total  pressure  is 

limited  to  values  greater  than  P  ,  ,  where  P  ,  is  a  material  parameter.  In 

min  min 

addition,  total  tensile  pressure  is  limited  as  specific  internal  energy,  I, 
approaches  Im»  the  value  at  which  strength  effects  are  no  longer  important. 
In  this  case  the  total  tensile  pressure  is  multiplied  by 


2.5  Max  (I  -  Max  (I,  0.61  ),0)/I 
m  mm 

This  formulation  insures  that  tensile  pressure  values  which  have  passed  the 

P_,  test  are  not  modified  unless  I  exceeds  0.6  I  . 
m  m  m 

The  total  j  essure  is  calculated  using  the  Mie- Grunoisen  form 


P  =  P0  (1  -  I>/2)  +  Tpl  +  Pq 
where  ,  the  Gruneisen  parameter,  is  given  by 


r  -  r„p„/p 


subject  to  the  restriction  that  T  cannot  exceed  2  F  ,  P  is  given  bv 

oo 


P 

o 


1.013x10* 

Io 


min 


(I. 


Io)  dynes/cm1 


where  Io  is  ambient  (1  atmosphere)  specific  internal  energy.  Pq  then  insures 
that  material  pressure  will  initially  be  1  atmos  phere.  The  total  solid  or 
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liquid  pressure,  P,  is  limited  by  Pyj^  aud  an  Im  factor  as  discussed  previous¬ 
ly.  In  addition,  it  is  limited  to  very  smell  positive  values  if  p  is'less 
then  -0.25  (the  theory  being  that  it  is  essentially  rubble  unable  to  support 
tension  at  these  low  densities). 

In  addition  to  P,  a  gas  pressure,  P_,  is  calculated  if  I  exceeds  I  .  This 

u  zn 

pressure  is  calculated  using  a  gamma  law  formulation 
P0  =  (y-l>  t  U-Im> 

PG  is  compared  to  the  solid/liquid  pressure,  P,  and  P  is  reset  to  the  maximum 
of  the  two  values. 

The  shear  modulus,  0,  is  calculated  from 

G  =  (G  p  -  (G  -G  )  (l-e^'^G) )  11— I/I  ) 
m  mo  m 

if  I  <  I 

m 

and 

G  =  0  if  I  2  I 

m 

where  Go#  G^  and  pfi  are  material  parameters. 

The  yield  strength,  Y,  is  calculated  from 

Y  =  Y  * Y’Hin  (A  +  BP,  Y  ) 

r  l  max 

where  Y  is  a  rubble  factor  set  as  follows 
r 

Y  «  1  if  p  >  -0.2 

r 

Y  «  0  if  p  <  -0.25 

r 

Y  »  20  (p  +0.25)  if  -0.25  <  p  <  -0.2 
r 

The  factor  insures  +hat  rubble  cannot  sustain  any  yield  strength.  Y^  is  an 
internal  energy  factor  and  is  given  by 

if  1  i  ** 

Uha  if  1  * 
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This  is  a  thermal  softening  factor.  The  term  A  +  BP  provides  a  Mohr-Coulomb 

yield  strength  which  reaches  saturation  at  Y  =  Y 
•  °  max 

Temperature  is  calculated  for  edit  purposes  by  arbitrarily  setting  ambient 
internal  energy  to  1  x  10®  ergs/gm  and  computing  a  C^  value  to  return  288°  k 
at  ambient  energy. 

The  attached  Table  V  lists  materials  which  utilize  this  equation  of  state, 
the  HULL  names  for  the  materials  and  the  material  quantities  built  into  HULL 
for  the  materials. 


TABLE  V 

HIGH  PRESSURE  GEOLOGIC  PROPERTIES 


Dry 

Dry 

Dry 

Saturated 

Moenkopi 

Kaibab 

Coconino 

Coconino 

MATERIAL 

Sandstone 

Limestone 

Sandstone 

Sandstone 

Limestone 

HULL  NAME 

MOENK 

KAIBAB 

DRYCOC 

WETCOC 

LIMEST 

Po(gm/cc) 

2.25 

2.30 

2.03 

2.35 

2.58 

pr(gm/cc) 

2.68 

2.70 

2.68 

2.35 

2.70 

K^(dynes/cma) 

60.  E9 

100. E9 

60.  E9 

225. E9 

500. E9 

K  (dynes/cm1) 

450. E9 

600. E9 

450. E9 

225. E9 

600. E9 

K  ( dynes /cm2) 

700. E9 

800. E9 

700. E9 

580. E9 

800. E9 

Kjjj  (dynes /cm2) 
^Lock 

ft1 

4.656E12 

2.942E12 

4.656E12 

2.22E12 

2.25E12 

0.029208 

0.03541 

0.0488 

0 

0.196 

0.65 

0.2963 

0.75373 

0.74468 

0.33333 

0.25 

0.4 

0.25 

0.3 

0.4 

Mq 

0.0015 

1.0 

0.0015 

0.0015 

1.0 

P  ,  (dynes/cm*) 

pmin 

(kynes/cm* ) 

-0.667E8 

-0.15E9 

-0.667E8 

-0.667E8 

-0.15E9 

13.55E9 

22.86E9 

23.08E9 

0 

130.87E9 

P„, (dynes/cm*) 
Gotdynes/cm*) 

397.14E9 

195.18E9 

468.17E9 

334.31E9 

221.44E9 

50. E9 

150. E9 

50.  E9 

45.  E9 

150. E9 

G^Cdynes/cm1) 

120. E9 

150. E9 

120. E9 

110. E9 

150. E9 

E  (ergs/gm) 

37. E9 

10. E9 

37. E9 

25.  E9 

10.  E9 

Am(dynes/cma) 

5.E7 

0.2E9 

5.E7 

5.E7 

0.2E9 

B 

0.75 

1.0 

0.75 

0.75 

1.0 

(dynes/cm*) 

3.E9 

3.E9 

3.E9 

3.E9 

3.E9 

CaX 

0 

0.5 

0.5 

0.5 

0.5 

0.5 

Y-l 

0.4 

0.4 

0.4 

0.4 

0.4 

■[ambient. 

(ergs/gm) 

1.E8 

1.E8 

1.E8 

1.E8 

1.E8 
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SECTION  V 

PROGRAM  SELECTION  AND  CONTROL 


HULL  System  elements  are  selected  and  configured  by  the  SAIL  code  with 
instructions  from  the  input  data  stream  and  PLANK.  The  default  option  direc¬ 
tory  on  the  SAIL  library  file  is  used  to  simplify  this  procedure  by  remem¬ 
bering  the  primary  option  set  values.  The  HULL  system  default  options  are  all 
keyed  from  the  installation  parameter  INST  and  have  the  value  and  installation 
definition  in  the  following  table. 


Value  of 
INST 

1 

2 

3 

4 

5 

6 

7 

8 
9 


Installation  Definition 

Air  Force  Weapons  Laboratory  (Kirtland  AFB,  NM) 

Air  Force  Armament  Laboratory  (Eglin  AFB,  FI,) 

McDonnel  Douglas  Astronautics  (Huntington  Bch,  CA) 

Army  Ballistic  Missile  Defense  Agency  (Huntsville,  AL) 
Science  Applications  Incorporated  (Huntsville,  AL) 
Atomic  Weapons  Research  Establishment  (UK) 

Army  Ballistic  Research  Laboratory  (Aberdeen,  MD) 
Vought  (Dallas,  TX) 

OTI  (Shalimar,  FL) 


Once  the  installation  parameter  is  set,  the  options  which  choose  the  computer 
system  specific  parameters  are  selected  automatically. 

After  the  correct  options  have  been  assembled  to  tailor  the  HULL  programs 
for  the  installation  machine  and  operating  system,  we  can  further  control  SAIL 
processing  to  produce  the  specific  form  of  KEEL,  HULL,  or  PULL  needed  to  run  a 
given  problem.  The  dump  tape  (TAPE4)  produced  by  KEEL  is  the  major  control 
for  setting  up  programs  HULL  and  PULL.  Most  of  the  input  to  KEEL  sets  these 
control  parameters  through  a  TAPE4  record  called  the  ZBLOCK.  The  names  of  the 
ZBLOCK  parameters,  the  meaning  of  their  various  values,  and  the  defaults  are 
given  below. 
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ZBLOCK  Name 


Values 


Definition 


(Alternate  name)  (Default) 

PROBLEM  (PROB)  Any  number  of  the  form  Arbitrary  problem  identifier 

XXXX.YYYY 


AREF 

True.,  (False.) 

Y(j=l)  boundary  reflective  (3D) 

ATMOS 

Integer  1 

tropical  atmosphere 

2 

temperate  atmosphere 

3 

arctic  atmosphere 

4 

exponential  atmosphere 

(5) 

constant  pressure  and  density 

BREF 

•TRUE. ,  (.FALSE.) 

bottom  boundary  reflective 

BURN 

Integer  (0) 

Y( j=l) ,  2D; 

Z(k=l) ,  3D 

uo  burn  routine 

1 

burn  front  advanced  by  fluid 

2 

state  (2D  only) 

programmed  burn  (2D  only) 

CODE 

Integer  (1) 

normal  hydrodynamic  -  elastic/ 

2 

plastic  formulation 

include  interactive  dust 

COLD 

( .TRUE. ) ,  .FALSE. 

radiation  cooling  if  false 

CYCLE 

Real 

number  of  time  steps  from 

DIMEN 

(2),  3 

initial  conditions 

number  of  spatial  dimensions 

DT 

Real 

time  step  (seconds) 

ELC 

Real 

energy  last  check,  current 
total  mesh  energy 
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EOS  1 

2 

(6) 


10 

ETH 

Real 

EXPAND 

Real  number  <1(.10) 

FAIL 

Integer  (0) 

1 

FLUXER 

Integer  0 

(3) 

FREF 

.TRUE.,  (.FALSE.) 

GEOM 

1 

(2) 

HOB 

Real  (0) 

IMAX 

Real  (100) 

IQ 

Real  (IMAX-1) 

ISLAND 

(0) 

1 

2 

JMAX 

Real  (200) 

Doan-Nickel  Air  equation  of  ftate 
Constant  gamma  equation  of  state 

Multi-material  Mie-Gruneisen 
equation  of  state 

tabulated  air  equation  of  state 

theoretical  total  mesh  energy 

mesh  expansion  fraction  for 
rezone 

no  failure  model 

simple  stress  criteria  (STRAIN=0) 
simple  stress  and  strain  criteria 
(STRAIN  =  1) 
no  fluxing 

volume  fluxing  with  species  mass, 
volume,  and  energy  retainod 

Y( j^JMAX)  boundary  reflective 
(3D) 

cartesian  coordinates 

cylindrical  coordinates  (2D  only) 

weapon  burst  height  in  kilometers 

number  of  zones  in  x  coordinate 
direction 

number  of  active  x  coordinate 
zones 

no  effect 

internal  mesh  reflective 
boundary 

reflective  boundary  with  rigid 
body  translation  (not  implemented 
in  3D) 

number  of  zones  in  y  coordinate 
direction 
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JQ 

Real  (JMAX-1) 

number  of  active  y  coordinate 
zones 

KMAX 

Real  (1) 

number  of  zones  in  z  coordinate 
direction 

KQ 

Real  (KMAX-1) 

number  of  active  z  coordinate 
zones  (3D) 

LREF 

(.TRUE.),. FALSE. 

X(i=IMAX)  boundary  reflective 

MAGFLD 

(0) 

no  magnetic  field 

1 

dipole  field  and  MHD  formulation 

METHOD 

1 

(first  order  in  time) 

(2) 

(second  order  in  time) 

MLC 

Real 

mass  last  check,  current 
total  mass  in  the  mesh 

NH 

Integer 

number  of  variables  per  zone 

(5  +  3*  NM  +  NHIST; 

2D) 

calculated  by  KEEL 

(6  +  3*  NM  +  NHIST; 

3D) 

NHIC 

Integer  (set  by  PLANK) 

number  of  words  in  core  to 
retain  mesh  variables 

NHIST 

Integer  (set  by  PLANK) 

number  of  fluxed  material  de¬ 
pendent  histories  per  cell 
(stress  components,  plastic 
work,  etc.) 

NM 

Integer  (1) 

number  of  materials 

NVARST 

NOP 

Integer  (15,  2D; 

22,  3D) 

Integer  (0) 

number  of  variables  per  station 

number  of  trace  particles  and 
stations 

NPLPB 

Integer  (calculated 
by  PLANK) 

number  of  planes  per  record 
block  on  TAPE4  (3D) 

NPP 

Integer  (2,  2D; 

3,  3-D) 

number  of  parameters  per  particle 
(usually  coordinates) 
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NROWPB 

Integer  (calculated 
by  PLANK) 

NSTN 

Integer  (0) 

PTSTOP 

Real  (1.E20) 

RAD 

Integer  (0) 

1 

2 

RADLOSS 

Real 

REZONE 

(0) 

1 

2 

3 

4 

5 

6 
7 

RREF  .TRUE.,  (.FALSE.) 

STRAIN  Integer  (0) 

1 

STRESS  0 

(1) 

SUME  Real 

TERAD  Real 

TLC  Real 


number  of  rows  per  record  block 
on  TAPE4 

number  of  stations 

problem  stop  time 

no  radiation  routines 

equilibrium  radiation  diffusion 

non-equil ibrium  radiation 
diffusion 

total  energy  radit'^d 
no  rezone 

horizontal  and  vertical  -  shock 
triggered 

particle  triggered  -  horizontal 
and  vertical 

horizontal  shock  follower 

vertical  shock  follower 

vertical  particle  follower 

same  as  5  with  shocks  diffused 

continuous  (only  rezone  tor  3-D) 

X(i=IMAX)  boundary  reflective 

no  strain  components  saved 

strains  saved  (only  if  STRESS  >0) 

no  stress  (elastic/plastic 
formulation  not  included) 

stress  calculation  included 

total  energy  lost  by  cooling 

total  energy  lost  by  radiation 

time  of  last  check  on  energy  and 
mass 
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TREE 

.TRUE.,  (.FALSE.) 

Y ( j=JMAX) ,  2D;  Z(k=KMAX) ,  3D 
boundary  reflective 

TTIME 

Real 

total  problem  central  processor 
time  (seconds) 

TTSTOP 

Real  (100.) 

total  central  processor  time 
stop  (hours) 

UREZ 

Real  (10.) 

x  coordinate  velocity  rezone 
trigger  (cm/sec) 

vise 

Integer  (0) 

no  explicit  artificial  viscosity 

1 

explicit  artificial  viscosity 
included 

VREZ 

Real  (10.) 

y  coordinate  velocity  rezone 
trigger  (cm/sec) 

VOIDS 

Integer  (0) 

no  voids  allowed 

1 

implicit  voids  permitted 

2 

explicit  voids  whenever  material 
VOID  is  found 

WORK 

Integer  (0) 

no  plastic  work  hardening 

1 

plastic  work  hardening  included 

WREZ 

Real  (10.) 

z  coordinate  velocity  rezone 
trigger  (3D)  (cm/sec) 

XOB 

Real  (0.) 

x  coordinate  of  burst  (3D) 
(kilometers) 

YOB 

Real  (0.) 

y  coordinate  of  burst  (3D) 
(kilometers) 

YGND 

Real  (-1.E20) 

location  of  ground  for  rezone 
down  control.  BREF  is  set  to 
zero  when  mesh  bottom  reaches 
this  value 

The  ZBLOCK  variables  are  identified  on  TAPE4  by  their  names  in  alphabetic 
form  along  with  their  corresponding  values.  KEEL  will  default  ZBLOCK  varia¬ 
bles  to  the  values  indicated  above.  If  a  different  value  is  desired  it  can  be 
input  through  the  KEEL  or  HULL  input  data  stream. 
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SECTION  VI 
PROGRAM  KEEL 


Program  KEEL  defines  the  initial  ZBLOCK,  the  mesh  coordinates,  and  the 
material  property  specifications  (P,  u,  v,  w,  I,  M,  V)  on  a  cell  by  cell  basis 
for  the  entire  mesh.  Input  to  KEEL  is  free  format  with  variable  names  and 
values  delimited  by  blank,  comma,  and  equal  (=).  The  KEEL  input  must  start 
with  the  word  KEEL.  This  allows  program  KEEL  to  search  for  and  find  its  data 
independent  of  its  location  in  the  input  file.  The  KEEL  input  is  logically 
divided  into  three  parts:  initial  ZBLOCK  definition,  mesh  generation,  and 
material  generation.  The  logic  flow  of  KEEL  is  indicated  by  Figure  9. 

A.  ZBLOCK  DEFINITION 

The  ZBLOCK  is  literally  defaulted  for  every  variable  name.  The  only 
variables  that  must  be  set  are  to  define  the  number  of  materials  and  their 
names.  It  is  usual  to  associate  a  problem  number  to  identify  the  type  of 
problem  being  run,  a  user  numeric  identifier,  and  to  use  the  numerals  to  the 
right  of  the  decimal  point  i<  indicate  the  sequence  or  series  of  the  calcula¬ 
tion  when  running  parameter  studies.  If  the  BOW  program  library  is  being  used 
to  keep  track  of  problems  and  their  associated  tapes,  it  is  imperative  that 
all  problem  numbers  be  unique.  BULL  problems  can  be  assigned  problem  numbers 
from  00000.0001  to  99999.9999. 

After  the  problem  number  is  defined  it  is  necessary  to  determine  the 
number  of  zones  to  be  used  in  each  coordinate  direction.  The  HULL  code  is  not 
sensitive  to  moderate  variations  in  cell  size  from  zone-to-zone  in  any  given 
coordinate  direction  (10  percent).  However,  it  is  not  tolerant  of  large 
aspect  ratios  (greater  than  2.5-3).  The  number  of  zones  should  be  chosen  to 
provide  enough  resolution  to  properly  define  material  interfaces.  Choosing 
the  best  number  of  zones  in  each  coordinate  is  always  a  compromise  most 
heavily  weighted  by  economic  considerations  since  the  computer  time  required 
to  run  a  calculation  is  inversely  proportional  to  the  size  of  the  smallest 
zone  and  directly  proportional  to  the  total  number  of  zones  in  the  calcula¬ 
tion.  Therefore,  decreasing  the  zone  size  by  a  factor  of  2  in  a  3D  calcula- 
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Figure  9. 


Logic  Flow  in  Program  KEEL. 
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tion  would  result  in  a  factor  of  16  increase  in  the  run  time.  The  number  of 
zones  in  the  x,  y,  and  z  coordinate  directions  is  set  by  IMAX,  JMAX,  and  KMAX 
respectively.  Two-dimensional  calculations  will  not  use  KMAX  (the  z  coordi¬ 
nate). 

The  next  order  of  business  is  to  select  the  number  of  materials  and  their 
names.  Table  I  (Page  28)  lists  the  materials  currently  recognized  by  the  HULL 
system.  Any  material  listed  in  Table  I  may  be  input  to  KEEL  but  present  code 
structure  limits  their  total  number  in  a  given  calculation  to  ten  (i.e., 
NM^.10).  Each  material  defined  must  be  assigned  a  number  from  1  to  NM.  The 
order  is  arbitrary. 

If  deviatoric  stresses  are  to  be  included,  set  STRESS  =  1.  Otherwise  the 
calculation  will  assume  all  materials  are  pure  hydrodynamic  media.  If  it  is 
dosired  also  to  retain  the  strain  components  for  each  cell,  set  STRAIN  =  1. 
This  is  necessary  if  STRAIN  or  P/Y  criteria  are  to  be  used  for  failure. 
Similarly,  if  material  failure  is  to  be  explicitly  modelled  it  is  accessary  to 
set  FAIL  =  1.  This  has  effect  only  for  STRESS  =  1.  If  FAIL  =  1,  STRESS  =  1, 
and  STRAIN  =  0;  failure  will  be  determined  by  a  default  stress  criteria,* 
maximum  principal  2  2.5  yield  stress.  If  STRAIN  =  1,  both  the  stross  criteria 
and  a  principle  strain  2  -7  will  initiate  failure.  Those  are  default  values 
and  should  be  used  with  caution  because  of  our  current  poor  understanding  of 
failuro  initiation. 

Rezone  choice  may  be  made  by  indicating  the  number  from  0  to  7  for  the 
particular  rezone  selected.  REZONE  =  0  implies  no  rezone.  This  can  be  al¬ 
tered  later  as  HULL  is  running.  Only  REZONE  =  7  is  supported  in  three  dimen¬ 
sions.  This  rezone  (in  2D  or  3D)  can  in  general  cause  a  continuous  transla¬ 
tion  of  the  mesh  coordinates  in  any  direction.  It  is  currently  implemented  to 
translate  the  entire  mesh  at  a  continuous  velocity  in  the  vertical  direction 
(y  direction  in  2D;  z  direction  in  3D).  The  translation  velocity  is  set  by 
the  zone  which  contains  station  number  one.  Use  of  the  continuous  rezone  may 
reduce  the  numerical  method's  implicit  viscosity  to  an  unacceptable  level.  In 
any  case  it  is  best  to  run  most  dynamic  solid  calculations  with  explicit 
artificial  viscosity  by  setting  VISC  =  1. 
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If  it  is  desired  to  run  with  Lagrangian  tracer  particles,  the  variable  NOP 
should  be  set  to  a  value  greater  than  the  number  of  particles  to  be  generated. 
The  particle  routine  will  count  and  set  NOP  to  the  actual  number  inserted.  If 
stations  are  to  be  generated,  NS1N  should  be  set  to  a  number  greater  than  the 
number  to  be  inserted  (  see  Page  72  for  the  minimum  value  for  NSTN).  Since 
all  stations  are  also  particles,  NOP  must  be  greater  than  the  number  of 
particles  plus  the  number  of  stations.  If  a  physical  tape  is  to  be  used  to 
store  station  data  (and  neither  BOW  nor  the  permanent  files  are  being  used) 
then  it  is  necessary  to  set  STAPE  =  VSN,  where  VSN  is  the  volume  serial  number 
of  the  tape.  This  will  only  function  on  CDC  machinery  and  results  in  the 
automatic  issue  of  the  tape  request  by  KEEL.  Otherwise,  you  must  specify  the 
station  tape  TAPE9  by  means  of  control  cards. 

The  stability  factor  STABF  (fraction  of  the  Courant  condition  time  step) 

—  7 

is  currently  defaulted  to  .75.  The  default  cycle  one  time  step  is  1  x  10 
seconds.  If  the  problem  type  being  run  has  had  difficulties  starting  during 
previous  HULL  runs,  DT  should  be  reduced. 

For  most  circumstances  an  initial  time  step  DT  of  10-*®  seconds  should 
suffice.  If  STABF  has  to  be  set  to  much  less  than  .5  something  is  seriously 
wrong,  check  your  initial  conditions. 

With  the  current  version  of  HULL  (Version  113),  ISLAND  =  2,  BURN,  and 
REZONE  =  1,  ° ,  i,  4,  5,  6  have  not  been  implemented  in  the  three-dimensionrl 
HULL  code. 

If  the  problem  being  generated  has  a  mesh  region  that  is  inactive  (no 
velocity,  high  pressures,  or  cxp’osive  burn  region)  some  initial  problem 
computer  time  can  be  saved  by  setting  IQ,  JQ  and  KQ  (3D).  This  will  cause  the 
mesh  in  that  part  of  the  coordinate  beyond  X(IQ),  Y(JQ),  or  Z(KQ)  to  be  left 
out  of  the  calculation  until  the  activity  is  such  that  the  velocity  in  any  of 
the  active  region  limit  exceeds  UREZ,  VREZ,  or  WREZ  for  the  X,  Y,  and  Z 
coordinate  directions  respectively.  This  feature  can  not  be  employed  in 
conjunction  with  the  continuous  rezone  (if  REZONE  =  7).  As  the  activity 
reaches  the  limits  of  the  active  region  the  corresponding  activity  index  IQ, 
JQ  or  KQ  will  be  incremented  one  zone  per  cycle  until  it  reaches  the  mesh 
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limit  I  MAX-1,  J  MAX-1,  or  KMAX-1.  The  defaults  for  IQ,  JQ,  and  KQ  are  IMAX-1, 
J MAX-1,  and  KMAX-1. 

The  last  input  in  this  section  of  KEEL  is  the  header  title.  This  is  80 
columns  of  whatever  information  the  user  wishes  to  use  to  identify  the  calcu¬ 
lation.  This  title  will  be  printed  at  the  beginning  of  the  output  for  each 
HULL  run  and  will  appear  in  the  titles  of  all  plots  produced  by  PULL.  The  80 
columns  of  information  must  be  preceded  by  a  card  image  with  the  word  HEADER. 

The  end  of  ZBLOCK  data  definition  is  signalled  by  the  word  MESH. 

An  example  of  this  part  of  the  KEEL  input  for  a  three-material  problem  to 
calculate  the  stress  induced  in  an  RHA  plate  impacted  by  copper  and  ignoring 
failure  could  be  set  up  by: 

KEEL  PR0B=1001 .001 

NM=3  AIR=1  RHA=2  CU=3  DIMEN=2 

STRESS=1  VISOl  DT=1.E10 

REZONE=0  I MAX-80  JMAX=120 

HEADER 

TEST  CALCULATION  TO  ILLUSTRATE  KEEL  INPUT 
MESH 

If  it  is  desired  to  insert  stations  to  record  str  s  rime  histories  at  10 
points  then 

N0P=11 ,  NSTN=11 

could  be  placed  between  or  on  any  card  preceding  HEADER.  The  almost  explicit 
(with  aid  of  the  ZBLOCK)  definition  makes  this  portion  of  KEEL  input  relative¬ 
ly  easy  to  use. 

B.  MESH  GENERATION 

The  next  logical  element  of  the  KEEL  input  is  definition  of  the  mesh 
coordinates.  This  part  of  the  KEEL  input  begins  with  the  keyword  MESH. 
Several  options  are  possible.  * ;c  most  simple  is  a  linear  mesh  in  all  coordi- 
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nate  directions.  This  is  set  up  by  establishing  the  mesh  boundaries  in  each 
coordinate  direction.  The  variable  names  are: 

10  coordinate  of  X(I=0) 

Y0  coordinate  of  Y(J=0) 

Z0  coordinate  of  Z(K=0)  (3D  only) 

XMAX  coordinate  of  X(I=IMAX) 

YMAX  coordinate  of  Y(J=JMAX) 

ZMAX  coordinate  of  Z(X=EHAX)  (3D  only) 

The  mesh  coordinates  are  established  for  the  boundary  of  each  zone  having 
indices  (I,  J,  X)  at  the  position  midway  between  its  center  and  zones  (1+1,  J, 
K),  (I,  J+l,  X),  and  (I,  J,  X+l).  Thus,  the  definition  of  X0,  Y0,  and  Z0  is 
to  locate  left,  back  and  bottom  of  the  first  cell.  Each  of  these  variablos 
must  be  paired  with  a  numeric  value.  There  are  restrictions.  In  two  dimen¬ 
sions  only  X0=0  is  allowed.  In  all  cases  X0,  Y0  and  Z0  are  defaulted  to  zero. 
The  zones  will  bo  of  constant  size  for  this  mesh  specification  with  zone  size: 

DX=(XMAX-X0) /IMAX 
DY=(YMAX-Y0) /JMAX 
DZ=(ZMAX-Z0) /ZMAX 

Thus  for  the  values  of  IMAX=80  and  JMAX=120  in  the  proceeding  sections  the 
input 


MESH  XMAX=40  Y0=30  YMAX=90 
would  produce  a  mesh  with  DX=DY=.5  centimeters. 

There  are  usually  economic  constraints  to  running  calculations  with  zones 
small  enough  for  sufficient  resolution  of  all  effects  everywhere  in  the  mesh. 
The  best  compromise  frequently  is  to  concentrate  a  fine  constant  grid  in  the 
region  of  most  interest  and  increase  the  zones  in  all  directions  away  from 
that  region.  This  is  especially  important  for  reasonable  treatment  of  solids 
since  perfectly  transmissive  mesh  boundaries  for  the  HULL  code  have  not  yet 
been  formulated.  Thus  if  the  calculational  mesh  is  too  small,  disturbances 
will  partially  reflect  from  its  boundary  and  propagate  to  the  interior  of  the 
mesh.  These  reflected  disturbances  may  disrupt  the  validity  of  the  calcula- 
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tion  in  the  region  of  interest.  This  mesh  generation  option  is  obtained  with 
the  keywords: 

MESH 


CONSTANT  SUBGRID  NX=NN 

X0=X 

XMAX=XX 

NY=MM 

Y0=Y 

YMAX=YY 

NZ=LL 

Z0=Z 

ZMAX=ZZ  (3D) 

NX,  NY,  and  NZ  are  the  number  of  zones  in  the  constant  subgrid  region  with 
boundaries  XO,  XMAX;  YO,  YMAX;  and  ZO,  ZMAX  respectively.  The  zone  dimensions 
within  this  region  are  defined  by 

DX=(XMAX-X0) /NX 
DY=(YMAX-YO)/NY 
DZ= ( ZMAX-ZO ) / NZ  (3D). 

NX,  NY  and  NZ  are  defaulted  to  IMAX,  JMAX,  and  KMAX  respectively.  The  re¬ 
maining  IMAX-NX,  JMAX-NY,  and  KMAX-NZ  zones  outside  the  constant  region  will 
be  increased  in  size  by  a  constant  factor  as  the  zone  index  increases  beyond 
NX,  NY.  or  NZ.  The  factors  are  RXPOS,  RXNEG;  RYPOS,  RYNEG;  and  RZPOS,  RZNEG 
for  zones  with  indices  greater  or  less  than  the  constant  subgrid  boundary  in 
the  X,  Y,  and  Z  coordinate  directious  respectively.  The  default  value  is  1,1 
for  all  these  factors,  i.e.,  a  10  percent  growth  rate  will  be  applied  for 
adjoining  zones  as  they  get  progressively  further  from  the  constant  subgrid 
region.  This  is  a  good  value  to  use  in  most  situations.  If  it  is  desired  to 
have  the  exterior  limits  of  the  mesh  stop  near  some  limit,  then  the  variables: 

X0LIM=X  XMLIM=XX 
Y0LIM=Y  YMLIM=YY 
Z0LIM=Z  ZMLIM=ZZ 

can  be  specified  along  with  a  large  value  of  the  corresponding  multiplying 
factors.  The  result  in  the  X  coordinate  direction  would  be  to  place  X(0)  2 
XOLIM  and  X(IMAX)  £  XMLIM.  The  other  coordinates  would  be  treated  in  a 
similar  manner  for  their  corresponding  variables. 


72 


If  the  above  lacks  flexibility  for  a  particular  application,  the  user  can 
specify  the  size  of  each  zone  with  the  variable  name/valne  pairs 

MESH  NX=il  DX=xl  NX=i2  DX  =  x2 

NY=jl  DY=yl  NY=j2  DY  =  y2 

NZ=kl  DZ=zl  NZ=k2  DZ  =  z2 

where  il,  i2  are  the  number  of  zones  with  zone  size  xl  for  the  first  i! 
zones,  and  zone  size  x2  for  the  next  i2  zones  until  the  sequence  sum  is  equal 
to  IMAX.  XO  defines  the  origin.  The  mesh  coordinates  Y  and  Z  will  be  handled 
in  an  analogous  way.  The  individual  value  of  each  occurrenceof  NX,  NY,  and  NZ 
can  be  from  one  to  IMAX,  JMAX,  and  KMAX  respectively.  Thus  for  a  brute  force 
mesh  setup  one  could  specify 


MESH  Y0=0 

X0=0 

NX=1 

DX=.l 

NX=2 

DX=.2 

NX=3 

DX=.3 

NY=2 

DY=.l 

NY=4 

DY=.5 

NY=2 

DY=1 

would  produce  mesh  coordinates  (assuming  IMAX=6  and  JMAX=8)  of: 

Xj=0,  .1,  .3,  .5,  .8,  1.1,  1.4;  i=integers  [0,  6] 

.1,  .2,  .7,  1.2,  1.7,  2.2,  3.2,  4.2;  j=integers  [0,  8], 

The  search  for  additional  mesh  data  will  terminate  w !  "h  the  keyword 
GENERATE. 

C.  MATERIAL  GENERATION 

The  KEEL  input  data  starting  with  the  word  GENERATE  is  used  for  inserting 
materials,  particles,  and  stations  in  their  desired  mesh  coordinates.  Al¬ 
though  the  nuclear  options  for  isothermal  spheres  and  other  constructs  still 
remain  in  the  code,  they  have  not  been  used  in  several  years  in  this  version 
of  the  code,  and  probably  do  not  execute  properly.  Their  use  is  covered  in 
previous  documents  (REF  4)  and  will  not  be  repeated  here.  Their  functions  can 
be  replicated  by  the  syntax  which  follows. 
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Material  can  be  inserted  into  the  mesh  by  KEEL  from  two  basic  sources — 
from  a  previous  HULL  or  EPIC3  calculation,  or  from  the  KEEL  input  data  stream. 
To  input  data  from  a  previous  calculation  the  keyword  FIRPJ.N  must  be  used. 
This  is  immediately  followed  by  the  word  HULL  or  EPIC3  depending  on  the  data 
source.  The  next  variable  name/value  pair  to  occur  is 

TAPE=VSN  or  PFN=permanent  file  name 

ID=permanont  file  identifier 


where  VSN  is  the  volume  serial  number  of  the  tape  which  contains  the  HULL  data 
and  PFN  and  ID  are  the  permanent  file  name  and  file  identifier  (ID),  respec¬ 
tively,  of  the  EPIC3  data.  HULL  data  is  supplied  by  any  previously  run  HULL 
calculation  in  a  standard  TAPE4  format.  EPIC3  data  is  provided  by  a  special 
dump  routine  which  is  described  in  Reference  30.  The  KEEL  input  at  this  point 
would  appear  as 


GENERATE 

FIREIN 

HULL 

TAPE=VSN 

FIREIN 

EPIC3 

PFN=pfn  ID=id. 

The  construct  for  defining  the  material  and  location  destination  in  the 
mesh  is: 


PACKAGE  material  geom  DELETE  geom  1  geom  2 

where  material  is  the  HULL  name  of  one  of  the  materials  defined  in  the  first 
part  of  the  KEEL  input  data  stream.  This  material  must  also  have  been  present 
in  the  calculation  which  produced  the  HULL  or  EPIC3  data  tape  being  input  or 
this  package  will  have  no  effect.  Geom  is  the  geometric  figure  and  its 
specifications  which  defines  the  region  in  the  mesh  that  is  to  be  filled  from 
the  input  data  tape.  The  geometi  c  figures  which  follow  (geom  1,  geom  2  to  a 
total  of  9)  are  deleted  from  the  first  geometric  configuration.  The  word 
DELETE  is  not  necessary,  since  the  appearance  of  a  second  and  subsequent 
geometry  specification  within  a  PACKAGE  implies  deletion  of  that  region  encom¬ 
passed  by  its  boundaries. 
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Each  of  the  geometric  figures  is  accompanied  by  a  list  of  parameters  which 
define  their  descriptions  and  position  in  the  mesh.  The  following  constructs 
are  for  describing  two-dimensional  material  generation, 

RECTANGLE  Xl=xl  X2=r2  Yl=yl  Y2=y2 

A  rectangle  with  vertex  coordinates 
(xl,  yl),  (xl,  y2),  (x2,  yl) ,  (x2,  y2). 

TRIANGLE  Xl=xl  X2=x2  X3=x3 

Yl=yl  Y2=y2  Y3=y3 

A  triangle  with  vertex  coordinates 
(xl,yl),  (x2,y2),  (x3,y3). 

CIRCLE  XC=xl  YC=yl  RADIDS=r 

a  circle  of  radius  r  with  center  at  (xl,yl). 

ELLIPSE  A=yl  B=b  C=xl  D=d 

An  ellipse  with  center  (xl,yl)  and  semiaxis  length  b 
in  the  y  coordinate  and  d  in  the  x  coordinate 

PARABOLA  A=a  B=b  C=c 

All  (xj.yj);  yj  1  a+b  (xj-c)3. 

HYPERBOLA  A=a  B=b  C=c  D=d 

All  (x^yj);  yA  2  a+b-v/l+(xi~c)"*/d* 

GNRLFIT  A=a  3=b  C=c  D=d  E=e  F=f 

All  (xpyj);  y^  2  ax^*  +  bx^,4  +  cxj*  +  dx^1  +  ex^  +  f. 

CURVE  TABLE  NPT  =  n  (xl,yl)  (x2,y2) - (xn,yn) 

For  n  100  points  and  xl  <  x2  <  x3  < - <xn 

All  (x^yj);  xm  2  x£  2  x^  and  y*  2  ya+<*r*a>t:: - - > 

^+1  xm 

The  three-dimensional  geometry  definitions  are  described  in  a  similar  way. 

Box  Xl=xl  X2=x2  Yl=yl  Y2=y2  Zl=zl  Z2=z2 

The  inside  of  the  prism  with  corners  (xl,yl,zl), 

(xl,yl,z2),  (tl,y2,zl),  (xl,y2,z2),  (x2,yl,zl), 

(x2, yl , z2) ,  ( t2,  y2 , zl) ,  <x2,y2,z2). 
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WEDGE  Xl=xl  Yl=yl  X2=x2  Y2=y2  Y3=y3  X3=x3  Zl=zl  Z2=z2 


The  inside  of  the  triangular  prism  with  base  in  the  XY 
plane  and  corners  (xl,yl),  (x2,y2),  (x3,y3)  in  the  planes 
Z=zl  and  Z=z2. 

PYRAMID  Xl=xl  Yl=yl  Zl=zl  X2=x2  Y2=y2  Z2=z2 

X3=x3  Y3=y3  Z3=z3  X4=x4  Y4=y4  Z4=z4 

The  inside  of  the  tetrahedron  defined  by  the  triples 
(xi,yi,zi);  i =1,2, 3, 4. 

CYLINDER  XC=xl  YC=yl  ZB=zl  ZT=z2  RADIO S=r 

The  interior  of  the  right  circular  cylinder  of 
radius  r  with  axis  parallel  to  the  Z  axis  and  centered 
in  the  XY  plane  at  (xl,yl).  The  cylinder  ends  are  at 
Z=zl  and  Z=z2 . 


ELLIPCYL  A=a  B=b  C=c  D=d  ZB=zl  ZT=z2 

The  interior  of  the  elliptical  cylinder  described  by 
the  parameters  a,  b,  c,  d  as  in  the  two-dimensional 
ellipse  with  ends  at  Z=zl,  Z=z2. 

PARCYL  A=a  B=b  C=c  ZB=zl  ZT=z2 

The  interior  of  the  parabolic  cylinder  analogous  to 
the  ellipse  above. 

HYPERCYL  A=a  B=b  C=c  D=d  ZB=zl  ZT=z2 

The  interior  of  the  hyperbolic  cylinder  analogous  to 
the  ellipse  above. 

SPHERE  XC=xl  YC=yl  ZC=zl  RADIUS=r 

The  interior  of  the  sphere  of  radius  r  at  (xl,yl,zl). 

RECTAROT  Xl=xl  X2=x2  Zl=zl  Z2=z2 

The  region  swept  by  the  rectangle  specified  rotated 
about  the  z  aiis. 

TRIROT  Xl=xl  Zl=zl  X2=x2  Z2=z2  X3=x3  Z3=z3 

The  region  swept  by  the  specified  triangle  rotated 
about  the  Z  axis. 

CIRCLEROT  XC=xl  ZC=zl  RADIUS=r 

The  region  swept  by  the  specified  cirtle  rotated  about 
the  Z  axis. 


ELLIPSEROT  AZ=a  BZ=b  CZ-c  DZ=d 


The  region  swept  by  the  specified  ellipse  rotated 
about  the  Z  axis. 

PARABLROT  AZ=a  BZ=b  CZ=c 

The  region  swept  by  the  specified  parabola  rotated 
about  the  Z  axis. 

GNRLFIT  A=a  B=b  C=c  D=d  E=e  F=f 

The  region  swept  by  the  polynomial  analogous  to  the 
two-dimensional  general  fit  rotated  about  the  Z  axis. 

CURVE  TABLE  NPT=N  xi ,  zl,  x2,  z2 - xn,  zn 

For  n£100,  the  region  swept  by  the  analog  of  the 
two-dimensional  curve  table  rotated  about  the  Z  axis. 

All  of  these  geometric  figure  specification  parameter  names  have  synonyms. 
The  following  names  in  each  set  are  equivalent. 


SET  1 

XO,  XI, 

XL, 

XLFFT,  XC, 

XCKNT,  XCNTR,  XCENTER 

SET  2 

YO,  Yl, 

YL, 

YB, 

YBOT, 

YBOTTOM,  YC, 

YCENT,  YCNTR,  Y CENTER 

SET  3 

ZO,  Zl, 

ZL, 

ZB, 

ZBOT, 

ZBOTTOM,  ZC, 

ZCENT,  ZCNTR,  Z CENTER 

SET  4 

X2,  XR, 

XRIGHT 

SET  5 

Y2,  YR, 

YF, 

YT, 

YTOP 

SET  6 

Z2,  ZR, 

ZT, 

ZTOP. 

An  additional  degree  of  flexibility  is  provided  by  allowing  each  of  the 
above  figures  to  be  displaced  or  rotated  about  each  of  the  coordinate  axis. 
The  displacement  distance  is  specified  by 

XCC=xl  YCC=yl  ZCOzl  (3D) 

where  the  occurrence  of  one,  two,  or  all  of  the  above  would  displace  the 
geometric  frame  of  reference  by  xl  in  the  X  coordinate  direction,  yl  in  the  Y 
coordinate  direction,  or  zl  in  the  Z  coordinate  direction. 

Rotation  of  the  geometric  figure  can  be  accomplished  by 

ANGLA=a  DANA=da  NDA=1 

ANGLB=b  DANB=db  NDB=m  (3-D  only) 

ANGLOc  DANC=dc  NDC=n  (3-d  only). 


Rotation  will  be  accomplished  with  respect  to  the  geometric  frame  of  reference 
defined  by  XCC,  YCC  and  ZCC.  A  positive  rotation  is  counter  clockwise  when 
viewed  along  the  axis  of  rotation.  All  of  these,  any  pair,  or  any  singlely 
defined  angle  could  be  defined  for  rotating  about  the  Z,  Y,  or  X  axis  by 
angles  a,  b  and  c  degrees  respectively.  The  DANA=da  NDA=1  portion,  when 
present,  indicates  a  replication  of  the  figure  1  times  at  increments  of  da 
degrees.  All  of  these  variables,  for  specification  of  displacement  or  rota¬ 
tion,  default  to  zero.  An  illustration  using  a  HULL  input  tape  is  shown 
below. 


GENERATE  FIREIN  HULL  TAPE=AA1 
PACKAGE  FE  RECTANGLE  XL=0  XR=10 

YB=0  YT=1 

YCC=- .5  XCC=-5  ANGLA=45 

would  result  in  producing  a  mask  (for  a  2-D  mesh)  in  a  slab,  of  thickness  1 
cm,  centered  on  the  origin,  at  an  angle  of  45  degrees  to  the  X  coordinate 
axis.  Data  would  bo  obtained  from  the  HULL  problem  on  TAPE  AA1  and  any  iron 
(FE)  properties  available,  as  seen  through  the  mask,  would  be  inserted  in  the 
new  mesh.  Subsequent  PACKAGE  input  geometries  would  obtain  data  from  the 
same  tape  until  terminated  by  ENDFIRE. 

The  most  common  usage  of  the  package  contruct  is: 

PACKAGE  mats* ial  [U=u,  V=v,  W=w,  I=i,  RHO=p] 
geometry  DELETE  geom  1  geom  2  -  geom  9. 

The  definition  of  the  material  and  geometry  parameters  is  as  before.  The 

parameters  enclosed  by  [ - ]  and  the  word  DELETE  are  optional.  Any  of  them 

may  appear.  They  are  intended  to  specify: 

u=X  component  of  material  velocity  (cm/sec) 
v=Y  component  of  material  velocity  (cm/sec) 
w=Z  component  of  material  velocity  (cm/ sec) 
i=internal  specific  energy  (ergs/gm) 
p=material  density  (gm/cml). 
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The  velocity  components  default  to  zero.  The  energy  and  density  default 
to  ambient  conditions  for  standard  temperature  and  pressure  as  defined  by  the 
equation  of  state  for  that  material. 

Figure  10  is  a  twodimensional  example  with  a  plot  illustrating  the 
effects  of  the  rotation. 

A  final  use  of  the  geometry  specifications  is  the  generation  of  Lagranginn 
trace  particles.  This  construct  is: 

PARTICLES  Geometry  DELETE  geom  1,  geom  2,  -  geom  9. 

The  description  geometry  and  geom  1  -  geom  9  are  treated  as  before,  but  are 
used  to  insert  one  Lagrangian  particle  in  each  cell  whose  center  lies  within 
the  region  defined.  The  word  DELETE  is  optional.  The  particles  are  counted 
as  they  are  inserted  to  keep  the  ZBLOCK  variable  NOP  properly  set. 

The  final  TEEL  input  set  is  station  generation.  This  is  of  the  form 

STATION  XS=xl  x2  x3  -  xl 

YS=yl  y2  y3 - ym 

ZS=zl  z2  z3 - zn 

any  number  of  such  sets  can  appear.  They  have  the  effect  of  producing  data 
collection  points  at  each  location  specified  by  xi,  yj,  zk  in  the  list  above. 
This  is  an  implied  loop  that  is  executed  as 

LOOP  z=zl,  z2 ,  z3  -  zn 

LOOP  y=yl,  y2,  y3  -  ym 

LOOP  x=xl,  x2,  x3  -  xl 

ENDLOOP  DEFINE  station  with  coordinates  (x,y,z). 

The  maximum  of  1,  m,  n  must  be  less  than  NSTN  as  input  to  KEEL.  The  use  of 
XS,  YS,  ZS  (3-D)  indicates  stations  that  are  to  remain  fixed  with  respect  to 
the  Eulerian  reference.  Lagrangian  stations,  stations  that  follow  the  flow, 
can  be  inserted  by  replacing  the  S  in  XS,  YS,  and  ZS  by  an  L.  A  station  can 
have  Lagrangian  or  Eulerian  behavior  in  each  of  the  coordinate  directions, 
i.e.,  XL,  YS,  ZL  is  permissible. 


7  9 


GENERATE 

PACKAGE 

FE 

CIRCLE  XC= 
RECTANGLE 

24  YC=i.2  RADIUS-- 1 . 2 
YB=1 . 2 

PACKAGE 

FE 

RECTANGLE 

XL=22 . 0  XR-25.2 

YB=  1.2  YT~40 

PACKAGE 

RHA 

V=1 . 5E5 

RECTANGLE  YB=-8, 4  YT--=0 

XCC=24  YCC=-2 . 6 

ANGLA=60 

Figure  10.  Example  GENERATE  Input  for  KEEL  and 
Resulting  Geometric  Configuration. 
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After  all  particles,  stations  and  packages  have  been  processed,  XEEL 
writes  out  the  mesh  data  on  TAPE4  and  the  station  data  on  TAPE9.  The  formats 
are  shown  in  the  accompanying  Figures  11  and  12.  KEEL  also  performs  chocks 
for  empty  zones  and  indicates  their  location  and  number.  A  column  xist  and 
material  and  particle/station  print  plot  is  also  produced  to  allow  the  user  to 
check  the  results  of  generation. 
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TAFS4  FORMAT 


(DUMP  TAPE) 


f  555.,  PROB,  CYCLE,  T  (Header  Synchronization) 

ZBLK(100,2)  (ZBLOCK  Description) 

X(IMAX) ,  iO,  Y(JMAX) ,  ZO,  Z(KMAX)  (Coordintte  Description) 

H  (NVARPB)l  (Cell  by  Coll  State  Variables) 

DDMP 

NBLKS 

J 

H(NOP*NPP)  IF  NOP>0  and  NOP*NPP  i  NVARPB  (Particle  and  Station 

Coordinates) 

LH(INT(NVARPB/NPP)*NPP)  N0P*NPP  >_  NVARPB 


Additional  dumps 


666. ,  666.,  666.,  666.  (End  of  tha  World  Trailer  ID) 

EOF 

EOF 


Figure  11.  HULL  Mesh  Dump  Tape  Format 
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TAPES*  FORMAT 


(STATION  TAPE) 


[STAT,  PROB 

Header  ZBLK  (100,2) 

Informa¬ 
tion  H(N0P*NPP)  N0P*NPP.£1020 

otherwise 

LH(NPP*INT(1G20/NPP))  for  number  of  records  needed 
fTime  (last  15  bits  holds  number  of  words  at  this  time) 


HYDRO  STRESS (2-D) 

P[STA-last  15  Density[STA-last  15  bits] 

bits]  X[Material  Code] 

Densi ty [Materia] 

Code  3 


STRESS  (3-D) 
Density[STA-last 
15  bits] 
XtMaterial  Code] 


u  u  u 

V  V  V 

w 

Repeated  W  6  ft 

for  next 

active  $  ^ 

station 

till  1024  ft 

word  SRR-P  SXX-P 

buffer  is  SRZ  SYY-P 

filled  SZZ-P  SZZ-P 

SH00P-P  SXY 

ERRL  SXZ 

EZZ  SYZ 

EHOOP  EXX 

Y  EYY 

I(NVARST>15)  EZZ 

EXY 

EXZ 

luYZ 

Y 

L  Z 


MATERIAL  CODE  -  bit  number  indicates  material  present  i.e.,  last  15  bits  for 
indicating  materials  bit  on  indicates  presence  of  th? t  ma¬ 
terial  =  material  3  and  5  is  00000  00000  10100. 

Figure  12.  HULL  Station  Tape  Format 
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SECTION  VII 
PROGRAM  HULL 


HULL  solves  the  finite  difference  equations  on  a  mesh  defined  by  KEEL  or 
a  previous  HULL  run.  The  mandatory  input  is  the  initial  conditions  contained 
on  TAPE4.  If  station  data  is  being  retained  (NSTN>0)  then  a  TAPE9  will  also 
be  input.  All  other  input  data  for  running  HULL  will  be  obtained  from  the 
input  data  stream.  The  HULL  input  mus  start  with  the  keyword  HULL.  All  data 
on  the  HULL  input  record  is  in  free  format,  delimited  by  blank,  comma,  and 
equal  (=).  The  HULL  input  record  is  logically  divided  into  three  parts: 
restart  information,  ZBLOCK  and  program  control,  and  solid  material  specifica¬ 
tions.  Figures  13  and  14  are  schematic  representations  of  the  data  flow  in 
program  HULL. 

A.  RESTART  DATA 

The  first  portion  of  the  input  record  immediately  following  the  keyword 
HULL  is  used  to  assure  the  correct  TAPE4  input  is  being  read.  The  parameter 
name/value  pair  PROB  =  XXXX.YYYY  is  compared  with  the  header  record  on  TAPE4. 
If  they  are  not  the  same  to  within  lO-*’  HULL  terminates  abnormally,  It  this 
test  is  passed,  HULL  checks  to  see  if  the  problem  time  and  cycle  number  on 
TAPE4  is  less  than  the  value  desired  for  restart.  If  they  arc  not,  the  TAPE4 
is  read  until  the  condition  is  met  or  the  end  of  file  as  designated  by  the 
HULL  trailer  is  reached.  If  end  of  file  is  reached,  HULL  backspaces  TAPE4  to 
the  beginning  of  that  dump  or  time  snap-shot  for  restart.  The  time  and  cycle 
for  restart  are  defaulted  to  l.xlO*0.  If  restart  at  a  particular  time  or 
cycle  is  desired,  the  parameter  name/value  pairs: 

T=t  CYCLE=c 

can  be  inserted  in  the  input  stream  immediately  after  the  PROB=XXXX.YYYY  pair. 
Only  T  or  CYCLE  need  be  input.  The  value  t  instructs  HULL  to  restart  on  the 
first  problem  dump  whose  problem  time  is  greater  t\an  or  equal  to  t.  Cycle 
starts  are  determined  in  the  same  way  except  the  comparison  is  doue  between 
the  value  c  and  CYCLE  as  it  appears  on  TAPE4.  To  be  used  for  restart,  all  of 
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MEMORY  RESIDENT  LOGICAL  DATA  MAPPING 


ROW/PLANE  N 
ROW/PLANE  N-l 
ROW/PLANE  N-2 
ROW/PLANE  N-3 


FIRST  WATCH 
SECOND  WATCH 
THIRD  WATCH 
FOURTH  WATCH 


PRINCIPAL  SUBROUTINE  LINKAGE 


FIRST  WATCH 


-EARTH 


-FIRE 

-WATER 


SECOND  WATCH 


THIRD  WATCH 


FOURTH  WATCH 


Figure  14.  Memory  Resident  Logical  Data  Mapping 

and  Principal  Subroutine  Lir.Kage  During 
In-Cycle  Processing. 
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the  names  PROB,  T,  and  CYCLE  must  appear  before  the  end  of  this  part  of  the 
HULL  input  record.  This  part  is  terminated  by  the  keyword  INPUT. 

A  typical  restart  segment  of  a  HULL  input  deck  is: 

HULL  PROB=1.001  CYCLE=105 
INPUT 

This  would  cause  HULL  to  read  the  TAPE4  header,  assuring  that  the  tape  con¬ 
tained  information  for  problem  number  1.001  a’  '  advance  TAPE4  to  the  first 
dump  whose  cycle  is  greater  than  or  equal  to  105.  If  cycle  100  is  the  last 
cycle  on  TAPE4,  it  would  be  positioned  to  restart  at  cycle  100. 

B.  Z BLOCK  AND  PROGRAM  CONTROL 

The  part  of  the  HULL  input  data  following  INPUT  and  preceding  the  end  of 
file  or  the  keyword  MATPR0P,  is  used  to  redefine  some  ZBLOCK  parameters  and  to 
specify  program  edit  control  and  tape  dumps.  Any  parameter  value  in  the 
ZBLOCK  can  be  modified  by  this  portion  of  the  input  stream  by  using  the 
appropriate  ZBLOCK  name/value  pair  in  free  format  as: 

NAME=value 

where  NAME  is  the  ZBLOCK  name  and  value  is  an  integer,  floating  point  or 
logical  Input  definition  as  defined  in  program  KEEL.  Care  should  be  exercised 
in  this  usage  since  it  is  possible  to  redefine  the  problem  number,  current 
time  or  cycle  by  inadvertantly  placing  the  keyword  PROB,  T,  or  CYCLE  after  the 
keyword  input  rather  thin  before.  Thus: 

HULL  PR0B=1 . 001 
CYCLE=24 
INPUT 
T=l.E-6 

would  result  in  starting  the  calculation  at  the  first  cycle  2  24  and  then 
redefining  the  problem  time  for  that  cycle  as  l.E-6  seconds.  The  following 
ZBLOCK  parameters  should  not  be  modified  since  they  will  change  the  process 
flow  of  the  code  or  the  meaning  of  an  edited  output  variable. 
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Dangerous  DULL  INPUT  ZBLOCK  Redefinitions 


SET  BY  KEEL  TO  SET  BY  HULL  TO 

DEFINE  THE  RUN  PROVIDE  INFORMATION 


CODE 

DIMEN 

IMAX 

JMAX 

KMAX 

NH 

NHIST 

NM 

NOP 

NROWPB 

NPLPB 

NPP 

NSTN 

NVARST 

PROB 


CYCLE 

ELC 

ETH 

MLC 

MTH 

SUME 

T 

TLC 

TTI.ME 

TTIME6 

TTIME7 


Resetting  the  values  of  IQ,  JQ,  KQ,  BURN,  FLUXER,  FAIL,  EOS,  or  STABF  will 
produce  results  which  may  be  harmful  to  the  health  of  the  calculation.  They 
should  also  be  used  with  care. 

There  are  several  variables  which  HULL  uses  for  control  during  a  calcula¬ 
tion  that  do  not  appear  in  the  ZBLOCK.  These  are  listed  along  with  their 
function  and  default  values  in  the  following  table: 

Non-ZBLOCK  HULL  Input  Control  Parameter* 


NAME 

FUNCTION 

DEFAULT  VAL 

DCYCST=n 

stop  calculation  in  n  cycles 

l.xlO20 

RTSTOP=t 

stop  calculation  after  this  run 
accumulates  t  CPU  hours 

l.xlO2 

DEBUG=L 

print  and  dump  every  cycle  if  L  is  . 

TRUE.  .FALSE. 

TIMES=n 

n=l  produces  36  dumps  per  decade  in 

n=2  produces  dumps  at  10ms,  30ms; 
every  .Is  from  ,1s  to  1,2s;  every 
.5s  from  1.5s  to  6.0s 

time  1 

n=3  produces  dumps  at  interval  DMPINT 
(commonly  used  for  most  conventional 
weapons  problems) 
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DMPINT=At 

VMIN=v 


ITRACE=i 

JTRACE=j 

KTRACE=k 


the  time  interval  At  for  producing  dumps  l,xlOao 

minimum  velocity  component  permitted;  1  cm/sec 

if  the  absolute  value  of  any  calculated 
velocity  component  is  loss  than  v  it  is 
zeroed 

cell  (i,j,k)  will  be  traced  by  printing  l.xlO4 
all  intermediate  values  when  SAIL  option 
DEBUG=2  is  set 


HEADER 

MRELER 


80  columns  from  the  following  card  will 
be  used  to  replace  the  ZBLOCK  header  title 

maximum  relative  error  allowed  per  cycle  l.xlO 
in  mass  or  energy 


This  section  of  HULL  input  terminates  with  the  keyword  MATPROP,  or  an  end-of- 
record. 


C.  MATERIAL  SPECIFICATIONS 

All  material  properties  for  HULL  are  defaulted  to  the  best  current  guess. 
If  a  calculation  involves  a  particular  material  specimen  for  which  experimen¬ 
tal  results  are  available,  it  is  in  general  best  to  replace  the  default  values 
of  the  yield  strength  parameters  by  their  experimentally  determined  values. 
The  parameters  that  can  be  changed  through  the  material  properties  specifica¬ 
tion  input  that  follow  MATPROP  are  in  a  material  dependent  format  that  can 
contain  the  following  optional  name/’vaiue  pairs. 

MATPROF 

MAT=n,  YLDST-yl,  YLDMAX=y2 ,  EPLAST=e,  IFRACl=ifl,  YFRACl=yfl, 

IFRAC2=if2,  YFRAC2=yf2 

ENDPROP 

If  the  keyword  MATPROP  appears  then  ENDPROP  must  appear.  The  specification  of 
each  set  of  specifications  for  a  given  material  must  be  preceded  by  MAT=n 
where  n  is  the  material  number  to  be  changed.  The  other  variable  name/value 
pairs  are  optional  but  names  must  be  spelled  correctly  or  the  routine  will 
search  forever  trying  to  recognize  them.  In  order  for  EPLAST=e  to  have  an 
effect,  the  ZBLOCK  must  contain  WORK  2.1>  otherwise  the  yield  surface  will  be 
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defined  by  YLDMAX.  When  WORK  is  on  (J>1),  plastic  work  (ep)  is  accumulated  for 
each  zone  which  is  occupied  by  a  solid.  The  yield  surface  is  then  defined  by 
Y  as  a  function  of  the  accumulated  plastic  strain  as  indicated  in  the  equation 
of  state  description.  Yield  is  always  a  function  of  internal  specific  energy. 
Thermal  softening  is  accomplished  by  a  trilinear  fit  as  previously  discussed 
in  the  equation  of  state  section.  The  thermal  factor  is  multiplied  times  the 
current  yield  surface  definition  (including  plastic  hardening)  to  obtain  a 
final  value  of  Y  (flow  stress)  against  wMch  the  yield  test  is  made.  If  you 
do  not  input  these  values,  the  default  values  from  Table  II,  page  30,  will  be 
used. 
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SECTION  VIII 
PROGRAM  PULL 


Program  PULL  accepts  the  calculational  mesh  data  from  TAPE4  or  the  sta¬ 
tion  data  from  TAPE9  and  produces  plotted  output  on  a  variety  of  graphics 
display  devices.  The  specific  plots  desired  cun  be  specified  by  the  user 
through  the  PULL  data  input  stream.  INPUT  to  PULL  is  free  format  with  input 
parameter  name/value  pairs  and  keywords  delimited  by  blank,  comma,  or  equal 
(=).  The  PULL  input  must  start  with  the  word  PULL.  Actually  two  different, 
program  sets  can  be  produced  by  SAIL  with  the  keyword  PULL.  The  first  is  the 
general  mesh  data  dump  plotter  which  produces  plots  from  the  HULL  TAPE4.  The 
second  program  set  is  for  sorting  and  plotting  the  station  tape  data  written 
by  HULL  on  TAPE9.  Program  PULL  data  flow  is  depicted  by  Figure  15. 

A.  DUMP  TAPE  PLOTS 

PULL  is  capable  of  making  three  basic  types  of  plots.  A  field  variable 
(density,  pressure,  energy,  etc)  versus  one  of  the  spatial  coordinates  (x,y, 
or  z)  is  called  a  histogram  and  is  abbreviated  for  use  in  the  keyword  struc¬ 
ture  as  HST.  Lines  of  constant  value  of  a  field  variable  in  a  two-dimensional 
representation  are  plotted  for  up  to  20  arbitrary  constants  or  contours. 
These  are  called  contour  plots  and  are  abbreviated  in  the  keyword  calls  by 
CONT.  Plots  of  vector  quantities  with  an  arrow  indicating  vector  direction 
and  the  arrow-length  indicating  vector  magnitude  can  be  plotted  for  all  vector 
quantities  such  as  velocity  and  gradients  of  scalars.  These  plots  are  avail¬ 
able  only  for  two-dimensional  planes.  Vector  plots  are  abbreviated  in  key¬ 
words  as  VECT.  When  gradients  of  scalars  are  selected,  both  a  vector  repre¬ 
sentation  of  the  gradient  and  contour  plots  of  the  gradient  magnitude  are 

produced.  Gradient  plots  are  abbreviated  in  keywords  as  GRAD. 

Histograms  can  be  obtained  in  any  coordinate  direction.  Those  in  the  x 
coordinate  direction  are  called  horizontal  histograms  and  are  abbreviated  by  H 
which  precedes  tht  designation  HST.  Histograms  in  the  y  or  z  coordinate 
direction  are  called  vertical  histograms.  They  are  abbreviated  by  V  for 
construction  of  keywords.  Thus  to  obtain  histogram  plots  of  any  field  vari- 
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READ  INPUT 
DECK 


Figure  15.  Program  PULL  Data  Flow. 
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able  XXX  the  keyword  would  be  XXXVHST  or  XXXHHST.  Contour  plots  are  called  by 
XXXCONT  and  gradient  plots  are  selected  by  XXXGRAD.  The  allowable  designa¬ 
tors,  XXX  and  their  meanings  are  tabulated  in  Table  VI. 


TABLE  VI 
SLOT  DESIGNATORS 


POSSIBLE  PLOTS 

XXX _ DEFINITION _ HHST  VEST  CONT  GRAD 


P 

D 

V 

W(3D) 

E 

D 

T 

M 

K 

RP 

RE 

RD 

SXX(3D) 

SYY(3D) 

SRR(2D) 

SZZ 

SHH ( 2D ) 

SRZ(2D) 

SXY ( 3D ) 

SXZ(3D) 

SYZ(3D) 

S2I 

SP1 

",P2 

EXX(3D) 

EYY(3D) 

ERS(2D) 

EZZ 

E0H(2D) 

ERZ(2D) 

EXY(3D) 

EXZ(3D) 

SYZ(3D) 

E2I 

EP1 

EP2 


HYDROSTATIC  PRESSURE  X 
X  COMPONENT  OF  VELOCITY  X 
Y  COMPONENT  OF  VELOCITY  X 
Z  COMPONENT  OF  VELOCITY  X 
SPECIFIC  INTERNAL  ENERGY  X 
DENSITY  X 
TEMPERATURE  X 
INDIVIDUAL  MATERIAL  DENSITY  X 
SPECIFIC  KINETIC  El*  ERGY  X 
PRESSURE  IN  ATMOSPHERES  X 


ENERGY  IN  UNITS  OF  AMBIENT  X 
DENSITY  RELATIVE  TO  AMBIENT  X 
X  COORDINATE  NORMAL  STRESS  X 
DEVIATOR 

Y  COORDINATE  NORMAL  STRESS 


DEVIATOR 

RADIAL  STRESS  DEVIATOR  X 
Z  OR  AXIAL  STRESS  DEVIATOR  X 
HOOP  DEVIATORIC  STRESS  X 
SHEAR  STRESS  X 
SHEAR  STRESS  X 
SHEAR  STRESS  X 
SHEAR  STRESS  X 
SECOND  STRESS  INVARIANT  X 
MAXIMUM  PRINCIPLE  STRESS  X 
MINIMUM  PRINCIPLE  STRESS  X 
STRAIN  FOR  SXX  X 
STRAIN  FOR  SYY  X 
STRAIN  FOR  SRR  X 
STRAIN  FOR  SZZ  X 
STRAIN  FOR  SHH  X 
STRAIN  FOR  SRZ  X 
STRAIN  FOR  SXY  X 
STRAIN  FOR  SXZ  X 
STRAIN  FOR  SYZ  X 
SECOND  STRAIN  INVARIANT  X 
MAXIMUM  PRINCIPAL  STRAIN  X 
MINIMUM  PRINCIPAL  STRAIN 


X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 


X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 


X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 


X 

X 

X 


X 

X 


X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 
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Thus  PCONT  would  be  the  allowable  keyword  for  selecting  pressure  contours. 

As  many  of  the  above  as  have  meaning  (i.e.,  SRR  is  not  available  for 
STRESS=0  run)  can  be  entered  for  a  single  PULL  run. 

Contour  plots  can  be  modified  by  a  list  of  variables  following  the  keyword 
to  specify  contour  values  to  be  plotted.  For  example: 

DC'INT  -11. 5  2345 

instructs  .PULL  to  draw  contour  lines  for  each  of  the  density  values  1,  1.5,  2, 
3,  4,  and  5  gm/cm3.  Contours  can  also  be  modified  by  ICTRS=n  and  CONYV=xl, 
x2,  ...,  xn  being  selected  for  the  contour  plot  requests  that  follow.  The 
value  of  n  can  be  from  2  to  20;  default  is  10.  The  values  of  the  x^  can  be 
any  real  number.  The  default  is  a  set  of  contours  equally  spaced  from  the 
rninimuui  to  the  maximum  in  n  steps.  In  addition,  contour  values  can  be  modi¬ 
fied  by  the  keywords  LOGV  and  LOGD.  LOGV  causes  all  contour  values  to  be  set 
to  values  of  lxlOn  and  3xl0n  from  the  maximum  value  down  for  all  decreasing 
integers  below  n  until  the  number  of  contours  is  exhausted.  LOGD  produces  the 
contours  of  Logarithm  base  10  of  the  data.  Contours  are  plotted  by  interpola¬ 
ting  the  data  and  smoothing  to  second  order.  This  can  be  turned  off  by 
including  the  keyword  NOCURV.  The  inclusion  of  the  keyword  COORD  will  cause 
the  zone  size  and  the  values  and  coordinates  of  the  minimum  and  maximum  to  be 
printed  on  each  plot.  If  particles  or  stations  were  included  in  the  calcul  - 
tion  (NOP^l),  the  particle  positions  can  be  presented  as  dots  on  the  contour 
plots.  This  is  obtained  by  including  the  keyword  PART  in  the  input  deck. 

Histograms  will  be  produced  at  the  mesh  boundary  lower  limit  as  a  default. 
By  specifying: 


PVI1ST  =  yl,  x2 ,  dx 

histograms  will  be  produced  along  all  vertical  columns  starting  at  xl  with 
increment  dx  to  x=x2. 

In  order  to  get  the  above  equivalent  plots  from  a  3-D  calculation  it  is 
necessary  to  specify  XPLANES,  YPLANES,  or  ZPLANES  followed  by  3  values  as 
histograms  were  described  above.  The  following  contour  specifications  will 
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then  be  satisfied  in  those  planes.  If  the  keyword  3DPL0TS  is  present  for  a  3D 
problem  then  the  nser  may  specify  RDCONT,  RECONT  or  RPCONT  to  plot  a  3D 
representation  of  the  selected  con'.our  values  when  viewed  in  the  plane  normal 
to  a  line  from  the  mesh  origin  to  the  position  EYE=x,  y,  z.  The  contours  will 
be  considered  opaque  for  this  representation. 

The  problem  times  at  which  plots  can  be  selected  are  set  by: 

CTIME=tl,  t2 ,  ...  tn 

DTIME=dt 

FTIME=tl 

LTIME=tn 

If  CTIME  is  present,  it  must  be  the  last  data  in  the  KEEL  input  stream.  Plots 
will  be  produced  for  all  times  tl,  t2,  ....  tn.  The  remaining  three  time 
s. lectors  will  produce  plots  from  all  dumps  after  tl  at  time  intervals  of  dt 
to  the  last  time  tn.  If  dt  is  not  specified,  all  dumps  between  tl  and  tn 
will  be  plotted.  If  the  keyword  STIME  is  included,  only  standard  times 
(TIMES=1  from  HULL  specification)  will  be  produced. 

U.  STATION  PLOTS 

Plots  of  the  station  parameters  can  be  obtained  by  following  the  keyword 
P'TLL  by  the  keyword  STATION.  Sorting  and  plotting  routines  are  written  to 
file  SAIL.  They  must  be  compiled  separately  and  in  sequence.  The  first 
compilation  must  be  executed  first  to  produce  an  unpacked  station  file  with 
all  time  data  for  station  one  first,  followed  by  all  data  for  station  two 
until  station  NSTN  is  reached.  If  the  HULL  calculation  that  produced  the  file 
was  a  pure  hydrodynamic  run  (STRESS=0),  then  the  execution  of  the  second 
compilation  will  plot  only  the  pressure,  density,  dynamic  pressure  and  impulse 
as  a  function  of  time.  The  ..  ,tions  to  be  plotted  can  be  determined  by 
FSTA=nl  LSTA=n2.  All  stations  from  nl  to  n2  will  be  plotted.  The  defaults 
are  nl=l  and  n2=NSTN.  If  STRESS21,  then  a  specific  set  of  plots  can  be 
requested.  These  are: 


95 


RADIAL 

AXIAL 

TOTAL 

STRESS 

STRESSDEV  - 

TENSOR 

FLOXM 


produce  plots  of  radial  motion 

produce  plots  of  axial  motion 

produce  plots  of  total  motion 

produce  plots  of  each  stress  component 

produce  plots  of  each  stress  deviator  component 

produce  plots  of  the  2nd  stress  invariant  and 
the  principal  stresses 

produce  a  plot  of  mass  flux 


We  have  thus  far  not  mentioned  the  use  of  the  problem  number  for  PULL 
input.  It  is  not  needed  unless  it  will  be  used  to  find  the  input  files  from 
the  tape  library  or  the  permanent  file  directory.  If  PROB=XXXXX.YYYY  is 
present  it  must  be  correct  for  this  run  or  PULL  will  terminate  abnormrlly. 


SECTION  EK 

SUNNING  A  CALCULATION 


The  preceding  sections  have  addressed  set  up  of  the  various  input  data 
streams  for  running  programs  in  the  HULL  system.  We  have  not  discussed  system 
considerations  because  of  the  frequent  changes  made  in  operating  systems. 
HULL  is  most  often  used  on  CDC  operating  systems.  We  therefore  follow  with  an 
example  control  card  and  data  deck  for  setting  up,  running,  and  plotting  a 
HULL  calculation.  We  will  assume  the  existence  of  HULLIB,  a  user  library 
containing  the  absolutes  of  programs  SAIL  and  PLANK,  and  relocatables  for  the 
HULL  prologue  routines  on  the  HULL  system  file. 


A.  KEEL  HON  DECK 


The  deck  below  will  set  up  50  (X  coordinate)  by  100  ty  coordinate)  zones 
with  DX=DY=.l  cm  from  XQ=0  and  Y0=~5,  respectively.  The  materials  will  be 
air,  iron  (FE)  and  copper  (CU).  Copper  will  be  placed  in  a  cylinder  of  radius 
1  cm  extending  from  Y=0  to  Y=5.  Iron  will  be  placed  in  a  slab  from  Y=.5  to 
Y=0.  Tna  iron  slab  is  given  an  initial  velocity  of  2.10s  cm/sec. 


ATTACH  (HULLIB, ID=(ocal) 

LIBRARY  (HULLIB) 

LABEL  (TAPE4 , W , T=0 , V SN=M INE ) 

PLANK. 

ATTACH  (OLD,HULLSYSTEMFILE, II>=local) 
SAIL. 

FTN  (I=SAIL, B=KEEL, 0PT=2 , L=0) 

KEEL. 

♦EOR  or  7/8/9  for  card  input  deck 
KEEL  PR0B=1 . 

IMAX=50  JMAX=100  DIMEN=2  STRESS=1 
NM=3  AIR=1  FE=2  CU=3 
REZ0NE=0  Header 

HULL  EXAMPLE  PROBLEM 
MESH  X0=0  XMAX=5  Y0=-5  YMAX=5 
GENERATE 

PACKAGE  AIR  RECTANGLE 

DELETE  RECTANGLE  XR=1  YB=0 
RECTANGLE  YB=.5  YB=0 
PACKAGE  CU  RECTANGLE  XR=1  YB=0 
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PACKAGE  FE  V=2.ES 

RECTANGLE  YB=-5  YT=0 
END 
*EOR 

SAIL  LINENO 

Change  cards  needed  by  this  run 
3.  HULL  RUN  DECK 


ATTACH  (HULLIB, ID=local) 

LIBRARY  (HULLIB) 

LABEL  (TAPE4 , R, VSN=MINE) 

PLANK. 

ATTACH  (OLD, HULLSYSTEMFILE, ID=local ) 
SAIL. 

FTN  ( I=SAIL, B=HULL, OPT=2 , L=0) 

HULL. 

*EOR 

SAIL  LINENO 
local  changes 
•EOR 

HULL  PROB=l .  CYCLE=0 
INPUT 

TIMES=3  DMPINT=5.E-6 
PTSTOP=50 . E-6 
MATPROP 

MAT=2  YLDMAX=15.E9 
MAT=3  YLDMAX=4.5E9 
ENDPROP 
*EOR 


C.  PULL  RUN  DECK 


ATTACH  (HULLIB, ID=local) 

LIBRARY  (HULLIB) 

LABEL  (TAPE4.R, VSN=MINE) 

PLANK. 

ATTACH  (OLD, HULLSYSTEMFILE, ID=local ) 

SAIL. 

FTN  (Is=SAIL,B=PULL, L=0,0PT=2) 

PULL. 

•EOR 

DCONT=7 . 8  7.9  8  8.1  8,2  8.3  8.6  8.7  8.8  9 

PVHST  WECT  COORD  NOCURV 

STIME 
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SECTION  X 

SENSE  SWITCH  CONTROL 


BOLL  can  be  controlled  for  set  np  by  SAIL  and  during  running  by  CDC  sense 
switch  settings.  If  the  sense  switches  are  set  before  PLANE  runs  they  have 
the  following  effects: 


SENSE  SWITCH  ON  BEFORE  PLANE  SDNS 
SENSE  SWITCH  NUMBER  EFFECT 

0  (All  switches  off)  ECS  store  of  the  mesh 

5  In-core  central  memory 

storage  of  the  mesh 

6  Disk  storage  of  the  mesh 

If  the  sense  switches  are  turned  on  after  HULL  or  PULI  have  started  the 
following  actions  will  be  taken: 

SENSE  SWITCH  ON  AFTER  HULL  OS  PULL  START 
SENSE  SWITCH  NUMBER  EFFECT 

1  Stop  at  end  of  this  cycle  (HULL) 
plot  this  dump  and  stop  (PULL) 

2  Route  current  output  to  printer  (both) 

3  Display  current  time,  cycle,  time  step 
and  estimate  of  running  time  on 
operators  console  (HULL) 

Display  time,  cycle,  and  type  of  plot 
being  made  (PULL) 

4  Change  dump  tapes  (HULL) 

Change  to  next  input  tape  (PULL) 

5  Start  double  buffering  of  a  disk  storage 
problem.  Increases  central  memory  but 
develops  I/O.  (HULL) 
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5  OFF  Stop  double  buffering,  return  unaeeded 

central  memory.  (HULL) 

6  Roll  ECS  out,  resume  by  operator  GO.  (HULL) 

The  sense  switches  can  be  toggled  on  by  the  operator  at  his  console  or 
through  the  use  of  control  cards  in  the  input  stream. 
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ATTN:  L.  K.  Goodwin 
Newport  Beach,  CA  92660 

3  General  Atomic  Company 

PO  Box  81608 
ATTN:  S.  M.  Snllivan 
F.  H.  Ho 
S.  Kwei 

San  Diego,  CA  92138 

1  Lockheed  Missiles  §  Space 

Company 

3251  Hanover  Street 
ATTN:  Org  5230,  Bldg.  201 
Mr.  R.  Robertson 
Palo  Alto,  CA  94304 

1  Lockheed  Missiles  and  Space 

Company 
PO  Box  504 

ATIN:  R.  L.  William* 

Dept.  81-11,  Bldg.  154 

Sunnyvale,  CA  94086 

1  Materials  Research  Laboratory, 

Inc. 

1  Science  Load 
Glenwood,  XL  60427 

2  McDonnel  1-Do eg .  Astronautics 

Company 

5301  Bolsa  Avenue 
ATTN:  Dr.  L.  B.  Greszczuk 
Dr.  J  Wall 

Huntington  Beach,  CA  92647 
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Copies  Organization  Copies  Organization 


5  Honeywell,  Ino, 

Government  and  Aerospace 
Prodncts  Division 
ATTN:  Mr.  J.  Blackburn 
Dr.  G.  Johnson 
Mr.  R.  Simpson 
Mr.  K.  H,  Doeringsfeld 
Dr.  D.  Vavrick 
600  Second  Street,  NE 
Hopkins,  MN  55343 

1  Hughes  Aircraft  Corporation 
ATTN:  Mr.  W.  Keppel 

MS  M-5,  Bidg.  808 
Tucson,  AZ  85706 

2  Kaman  Sciences  Corporation 
P.  0.  Box  7463 

ATTN:  Dr.  P.  Snow 

Dr.  D.  Williams 
Colorado  Springs,  CO  80933 

1  Pacific  Technical  Corporation 

460  Ward  Drive 
ATTN:  Dr.  F.  K.  Feldmann 
Santa  Barbara,  CA  93105 

1  Rockwell  International 

Missile  Systems  Division 
ATTN:  A.  R.  Glaser 
4300  E.  Fifth  Avenue 
Columbus,  OH  43216 

3  Schumberger  Well  Services 
Perforating  Center 

ATTN:  J.  E.  Brooks 
J.  Brookraan 
Dr.  C.  Aseltine 
PO  Box  A 

Rosharon,  TX  77543 


1  New  Mexico  Institute  of 

Mining  and  Technology 
ATTN:  TERA  Group 
Socorro,  NM  87801 

1  Northrup  Corporation 

3901  W.  Broadway 
ATTN:  R.  L.  Ramkuraar 
Hawthorne,  CA  90250 

1  Nuclear  Assurance  Corporation 
24  Executive  Park  West 

ATTN:  T.  C.  Thompson 
Atlanta,  GA  30245 

2  Orlando  Technology,  Inc. 

PO  Box  855 

ATTN:  Mr.  J.  Osborn 
Mr.  D.  Matuska 
Shalimar,  FL  32579 

1  US  Steel  Corporation 

Research  Center 
125  Jamison  Lane 
Monroeville,  PA  15146 

1  VPI  p  SU 

106C  Norris  Eall 
ATTN:  Dr.  M.  P.  Kamat 
Blacksburg,  VA  24061 

2  Vought  Corporation 
PO  Box  225907 
ATTN:  Dr.  G.  Hough 

Dr.  Paul  M.  Kenner 
Dallas,  TX  75265 

1  Westinghouse,  Inc. 

PO  Box  79 

ATIN:  J.  Y.  Fan 

W.  Mifflin,  PA  15122 
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No.  of 

Copies  Organization 

1  Science  Applications,  Inc. 

101  Continental  Boulevard 
Suite  310 

El  Segundo,  CA  90245 

1  Ship  Systems,  Inc. 

11750  Sorrento  Valley  Road 
ATTN:  Dr.  G.  G.  Erickson 
San  Diego,  CA  92121 

1  Systems,  Science  and  Software 
PO  Box  1620 

ATTN:  Dr.  R.  Sedgwick 
La  Jolla,  CA  92037 

2  TRW  Systems  Group 

One  Space  Park,  Rl/2120 
ATTN:  D.  Ausherman 
M.  Bronstein 

Redondo  Beach,  CA  90278 

1  United  Technologies 

Roscarch  Center 
438  Weir  Street 
ATTN:  P.  R.  Fitzpatrick 
Glastonbury,  CT  06033 

1  University  of  Arizona 

School  of  Engineering 
ATTN:  Dean  R.  Gallagher 
Tucson.  AZ  85721 

1  University  of  California 

504  Hilgard  Avenue 
ATTN:  Dr.  H.  Ziv 
Los  Angeles,  CA  90024 


No.  of 

Copies  Organization 


1  Drexel  University 
Department  of  Mechanical  Engr. 
ATTN:  Dr.  P.  C.  Chou 

32d  and  Cnestnut  Streets 
Philadelphia,  PA  19104 

3  Southwest  Research  Institute 
Dept,  of  Mechanical  Sciences 
ATTN:  Dr.  U.  Lindholm 

Dr.  W.  Baker 
Dr.  R.  Whito 
8500  Culebra  Road 
San  Antonio,  TX  78228 

4  SRI  International 
333  Ravenswood  Avenue 
ATTN:  Dr.  L.  Seaman 

Dr.  L.  Curran 
Dr.  D.  Shoe key 
Dr.  A.  L.  Florence 
Menlo  Park,  CA  94025 

2  University  of  Arizona 

Civil  Engineering  Department 
ATTN:  Dr.  D.  A.  DaDoppo 
Dr.  R.  Richard 
Tucson,  AZ  85721 

1  University  of  Oklahoma 

School  of  Aerospace, 

Mechanical  and  Nuclear 
Engineering 
ATTN:  Dr.  C.  W.  Bert 
Norman,  OK  73069 
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No.  of 

Copies  Organization 

1  University  of  California 
Department  of  Physics 
ATTN:  Dr.  Harold  Lewis 
Santa  Barbara,  CA  93106 

2  University  of  California 

College  of  Engineering 

ATTN:  Prof.  W.  Goldsmith 

Dr.  A.  G.  Evans 
Berkeley,  CA  94720 

2  University  of  Delaware 

Department  of  Mechanical 

Engineering 

ATTN:  Prof.  J.  Vinson 
Prof.  B.  Pipes 
Newark,  DE  19711 

1  University  of  Denver 

Denver  Research  Institute 
ATTN:  Mr.  R.  F.  Recht 
2390  S.  University  Blvd. 
Denver,  CO  80210 

2  University  of  Florida 

ATTN:  Dr.  R.  L.  Sierakowski 

Dr.  L.  E.  Malvern 

Gainesville,  FL  32061 


No.  of 

Aberdeen  Proving  Gronnd 

Dir,  USAMSAA 
ATIN:  DRXSY-D 

DRXSY-MP,  H.  Cohen 
Cdr,  USATECOM 
ATTN:  DRSTE-TO-F 
Dir,  USAAPG 

ATTN:  S.  Keithely,  »?rD 
Dir,  USACSL,  Bldg.  E3516,  EA 
ATTN:  DRDAR-CLB-PA 
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USER  EVALUATION  OF  REPORT 


Please  take  a  few  minutes  to  answer  the  questions  below;  tear  out 
this  sheet,  fold  as  indicated,  staple  or  tape  closed,  and  place 
in  the  mail.  Your  comments  will  provide  us  with  information  for 
improving  future  reports. 

1.  BRL  Report  Number _ 

2.  Does  this  report  satisfy  a  need?  (Comment  on  purpose,  related 
project,  or  other  area  of  interest  for  which  report  will  be  used.) 


3.  How,  specifically,  is  the  report  being  used?  (Information 
source,  design  data  or  procedure,  management  procedure,  source  of 
ideas,  etc.)  _ 


4.  Has  the  information  in  this  report  led  to  any  quantitative 
savings  as  far  as  man-hours/contract  dollars  saved,  operating  costs 
avoided,  efficiencies  achieved,  etc.?  If  so,  please  elaborate. 


5.  General  Comments  (Indicate  what  you  think  should  be  changed  to 
make  this  report  and  future  reports  of  this  type  more  responsive 
to  your  needs,  more  usable,  improve  readability,  etc.) _ 


6.  If  you  would  like  to  be  contacted  by  the  personnel  who  prepared 
this  report  to  raise  specific  questions  or  discuss  the  topic, 
please  fill  in  the  following  information. 

Name: 


Telephone  Number: 
Organization  Address: 


FOLD  HERE 


Director 

US  Army  Ballistic  Research  Laboratory 
Aberdeen  Proving  Ground,  MD  21005 


OFFICIAL  BUSINESS 

PENALTY  FOR  PRIVATE  USE,  S300 


BUSINESS  REPLY  MAIL 

FIRST  CLASS  PERMIT  NO  12062  WASHINGTON, DC 
POSTAGE  WILL  BE  PAID  BY  DEPARTMENT  OF  THE  ARMY 


Director 

US  Army  BallisticHlesearch  Laboratory 
ATTN:  DRDAR-TSB 

Aberdeen  Proving  Ground,  MD  21005 


NO  POSTAGE 
NECESSARY 
IF  MAILED 
IN  THE 

UNITED  STATES 


FOLD  HERE 


